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4.  INTRODUCTION 

The  goal  of  the  research  project,  “Time-Resolved  and  Spectroscopic  Three-Dimensional  Optical 
Breast  Tomography,”  is  to  develop  a  safe  and  affordable  breast  cancer  detection  modality  that 
uses  noninvasive  near-infrared  (NIR)  light  for  imaging,  locating  and  diagnosing  tumors  in 
human  breast  with  high  resolution  and  specificity.  The  research  focuses  on  the  development  and 
application  of:  (a)  continuous  wave  (CW)  and  time-resolved  experimental  approaches  for 
probing  the  target  using  NIR  light,  and  (b)  image  reconstruction  and  target  localization 
algorithms  for  obtaining  target  location  and  generating  three-dimensional  (3-D)  tomographic 
images. 

Significant  progress  has  been  made  in  both  the  thrust  areas  during  the  reporting  period  (May  1, 
2006  -  April  30,  2007)  covered  by  this  report. 

5.  BODY 

The  tasks  performed  and  the  progresses  made  during  the  current  reporting  period  are  as  follows: 

•  Use  of  optical  tomographic  imaging  using  independent  component  analysis  (OPTICA) 
developed  during  the  first  reporting  period  to  image  targets  in  ex  vivo  breast  tissue 
specimens, 

•  Extension  of  OPTICA  for  obtaining  cross  sectional  image  of  targets  in  turbid  media, 

•  Development  of  a  forward  model  for  scattered  light  intensity  distribution  in  the 
backscattering  geometry;  and 

•  Development  of  an  approach  for  determination  of  light  absorption  and  scattering 
characteristics  and  anisotropy  factor  of  a  highly  scattering  medium  using  backscattered 
circular  polarized  light. 

We  provide  a  brief  outline  of  our  accomplishments  in  these  areas,  and  refer  to  appended 
publications  for  detailed  description  where  applicable. 

5.1.  Imaging  and  localization  of  fluorescent  targets  embedded  in  ex  vivo  breast  tissue 
specimen 

12 

We  used  the  experimental  arrangements  ’  that  were  set  up  and  tested  [Technical  Objective 
(TO)  1,  Task  1 ]  during  the  first  reporting  period  and  detailed  in  the  first  report  to  carry  out 
imaging  and  localization  measurements  on  realistic  breast  models  assembled  using  ex  vivo  breast 
tissues  (TO  1,  Task  2  and  TO  2,  Task  4).  Fluorescent  contrast  agents  and  native  fluorescence 
have  the  potential  to  enhance  contrast  of  tumors  compared  to  the  normal  tissue.  So,  one  of  our 
aims  was  to  explore  the  feasibility  of  imaging  and  locating  a  fluorescent  target  inside  the  model 
breast. 

The  ex  vivo  female  breast  tissue  specimen  was  a  100  mm  X  100  mm  X  26  mm  slab  held 
inside  a  transparent  plastic  box.  The  fluorescent  target  was  a  glass  sphere  (outer  diameter  ~4  mm 
and  inner  diameter  3.2  mm)  filled  with  a  solution  of  indocyanine  green  (ICG)  dye  in  water.  The 
water  solution  of  ICG  had  an  absorption  coefficient  of  1 1.5  cm'1  at  784  nm,  and  fluoresced  over 
the  790-966  nm  spectral  range  with  a  peak  at  825  nm. 
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A  200- pm  optical  fiber  delivered  a  784-nm,  300  mW  continuous-wave  beam  from  a  diode 
laser  for  sample  illumination.  The  beam  was  collimated  to  a  1-inm  spot  onto  the  entrance  face 
(the  ‘source  plane’)  of  the  slab  sample.  Multiple  source  illumination  was  realized  in  practice  by 
step  scanning  the  slab  samples  across  the  laser  beam  in  an  1 8  X  15  x-y  array  of  grid  points  with  a 
step  size  of  2.5  mm  using  a  computer  controlled  translation  stage.  The  signal  from  the  opposite 
face  of  the  sample  (the  ‘detection  plane’)  was  passed  through  a  narrow-band  interference  filter 
centered  at  830  mn  (FWHM  ~  10  mn,  50%  transmission)  to  block  the  scattered  784-nm  pump 
light.  A  camera  lens  collected  the  fraction  of  filtered  fluorescence  around  830  mn  and  projected 
it  onto  the  sensing  element  of  a  cooled  16-bit,  1024  x  1024-pixel  charged  couple  device  (CCD) 
camera.  Each  illuminated  pixel  of  the  CCD  camera  could  be  regarded  as  a  detector. 

Each  raw  image  was  then  cropped  to  select  out  the  infonnation-rich  region,  and  binned  to 
enhance  the  signal-to-noise  ratio.  All  the  binned  images  corresponding  to  illumination  of  the  grid 
points  in  sequence  were  then  stacked,  and  used  as  input  for  independent  component  analysis.  The 
details  of  the  analysis  method,  theoretical  formalism,  target  localization  algorithm,  and 
experimental  arrangement  have  been  published,4  and  are  presented  in  Appendix  1 :  “ Three- 
dimensional  localization  and  cross  section  reconstruction  of  fluorescent  targets  in  ex  vivo  breast 
tissue  using  independent  component  analysis,"  Appl.  Phys.  Lett.  89,  133902  (2006). 

The  key  result  is  that  the  x-y-z  location  of  the  target  was  detennined  to  be  (30  mm,  27  mm, 
15  mm),  while  the  known  position  was  (30  mm,  27  mm,  16  mm).  The  lateral  position  of  the 
target  was  determined  exactly,  while  the  axial  position  was  determined  to  within  1  mm,  which  is 
quite  remarkable! 

5.2.  Extension  of  OPTICA  for  obtaining  cross  sectional  image  of  targets  in  turbid 
media 

Another  important  development  during  the  current  reporting  period  has  been  the  extension  of 
OPTICA  to  estimate  the  size  and  shape  of  the  target  (such  as,  tumor)  inside  the  turbid  medium 
(such  as,  breast).  OPTICA  was  initially  developed  ’  to  obtain  accurate  three-dimensional 
position  infonnation.  The  extension  of  the  formalism  to  generate  size  and  shape  information 
provides  inverse  image  reconstruction  ability  (TO  3)  and  its  application  to  targets  inside  realistic 
human  breast  models  ( T04 ,  Task  9)  is  a  crucial  step  towards  in  vivo  breast  imaging. 

In  this  extension  of  OPTICA,  the  size  and  shape  of  the  y'-th  target  within  the  sample  is 
estimated  from  a  back-projection  of  the  signal  due  to  the  target  from  the  detection  plane  onto  the 
target  plane.  Independent  component  analysis  (ICA)  provides  an  estimate  of  the  contribution 
Umj(rd,  i "j,  (o)  of  the  y-th  target  to  the  overall  signal  Um(  rt/,  iy,  co)  due  to  all  the  targets,  where  r y 
and  ly,  are  the  positions  of  the  detector  plane  and  the  y'-th  target  respectively.  The  mathematical 
formalism  for  obtaining  the  cross  sectional  image  is  summarized  in  Eq.  (4)  and  Eq.  (5)  of 
Reference  1  (Appendix  1),  which  also  provides  examples  of  the  application  of  the  formalism  to 
realistic  breast  model  (TO  4,  Task  9). 

The  efficacy  of  the  extended  OPTICA  formalism  was  tested  in  the  experiment  involving  ex 
vivo  human  breast  tissue  described  in  Section  5.1  above.  The  full  width  at  half  maximum  height 
(FWHM)  of  the  fluorescent  sphere  was  determined  to  be  4  nun  along  x  and  y  coordinate  axes 
which  agree  with  the  known  diameter  of  the  fluorescent  sphere  target. 
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5.  3.  Localization  and  cross  section  imaging  of  a  tumor  in  a  model  breast  sample 

The  OPTICA  formalism  along  with  its  extension  for  obtaining  cross  section  images  was  then 
applied  to  a  model  human  breast  sample  assembled  using  ex  vivo  breast  tissues  (TO I,  Task  2;  TO 
4,  Task  9).  In  this  application  the  difference  in  light  scattering  characteristics  of  the  tumor  and 
different  breast  tissue  components  was  used  for  target  detection.  The  details  of  the  experimental 
arrangement,  theoretical  formalism,  and  results  are  presented  in  Appendix  2:  “ Optical  high 
resolution  cross  section  imaging  of  a  human  breast  model  using  independent  component  analysis 
in  CW and  time  domain ,”  Proceedings  of  SPIE,  Vol.6434,  643416  (2007). 

The  sample  breast  had  a  dimension  of  42  mm  X  33  mm  X  30  mm.  A  tumor  of  dimension  5 
mm  X  5  mm  X  3  mm  was  placed  at  the  mid-plane  of  the  breast  sample.  The  sample  was 
illuminated  in  by  CW  784  nm  beam  from  a  diode  laser  delivered  by  a  200-pm  optical  fiber. 
Multiple  source  illumination  was  realized  by  step  scanning  the  sample  in  2-mm  steps  over  an 
array  of  22  X  16  grid  points  along  the  horizontal  (x)  and  vertical  (y)  directions  across  the  laser 
beam.  After  optical  measurements,  the  sample  was  transferred  to  our  collaborators  at  the  New 
York  Eye  and  Ear  Infirmary  for  pathological  study  and  correlation  ( T02 ,  Task  5). 

The  key  results  are  as  follows,  (a)  The  location  of  the  tumor  was  determined  to  within  1  mm 
in  all  three  dimensions  (details  in  Fig.  2,  Appendix  2);  (b)  the  FWHM  of  the  tumor  is  estimated  to 
be  ~  11  mm  (details  in  Fig.  3,  Appendix  2);  (c)  Two  other  structures  with  weaker  scattering 
coefficients  were  identified  in  the  mid-plane  of  the  sample.  Since  nonnal  component  of  the 
breast  tissue  specimen  was  mainly  adipose  tissue,  we  tentatively  identified  these  two  components 
as  the  fibro-glandular  tissues.  Comparison  with  the  pathology  results  further  confirmed  the 
identity  of  the  tumor  and  the  fibro-glandular  tissues. 

5.4.  Development  of  a  forward  model  for  scattered  light  intensity  distribution  in  the 
backscattering  geometry 

While  the  approaches  highlighted  above  use  experimental  data  in  transillumination 
geometry,  an  approach  based  on  backscattering  geometry  will  be  advantageous  for  targets  closer 
to  the  input  plane,  and  for  obtaining  depth  location.  We  developed  an  analytical  cumulant 
solution  of  both  the  time-dependent  scalar  radiative  transfer  equation  (SRTE),  and  the  vector 
radiative  transfer  equation  (VRTE)  in  an  infinite  uniform  medium.5’6  The  solution  of  SRTE  is 
suitable  for  describing  scattered  photon  density  distribution  for  an  arbitrary  phase  function  as 
long  as  it  is  a  function  of  the  scattering  angle.  The  numerical  results  of  light  propagation  from  the 
analytical  solution  are  in  a  good  agreement  with  that  of  Monte  Carlo  simulation.7 

However,  the  SRTE  cannot  account  for  the  polarization  of  the  incident  and  scattered  light 
beams.  Polarization  is  an  important  parameter  and  its  evolution  as  a  polarized  beam  of  light 
interacts  with  a  scattering  medium  with  embedded  targets  can  provide  useful  imaging 
information.  To  capitalize  on  the  advantages  that  polarization  has  to  offer,  we  used  the  analytical 
solution  of  the  VRTE  to  study  the  backscattering  of  circularly  polarized  light  pulses  from  an 
infinite  uniform  scattering  medium  as  a  function  of  the  helicity  of  the  incident  light  and  size  of 
scatterers  in  the  medium  (TO 3,  Task  6).  The  formalism  considers  a  polarized  short  pulse  of  light 
incident  on  the  scattering  medium,  uses  the  analytical  cumulant  solution  of  the  VRTE  with  the 
phase  matrix  obtained  from  the  Mie  theory  to  calculate  the  temporal  profile  of  scattered 
polarized  photons  for  any  position  and  any  angle  of  detection.  The  details  of  the  formalism, 
results,  and  comparison  with  experimental  results  are  presented  in  Appendix  3:  “Analytical 
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cumulant  solution  of  the  vector  radiative  transfer  equation  investigates  backs cattering  of 
circularly  polarized  light  from  turbid  media,  ”  Phys.  Rev.  E  74,  056605  (2006). 

The  results  show  that  circularly  polarized  light  of  the  same  helicity  dominates  the 
backscattered  signal  when  scatterer  size  is  larger  than  the  wavelength  of  light.  For  the  smaller 
particle  size,  light  of  opposite  helicity  makes  the  dominant  contribution  to  the  backscattered 
signal,  as  shown  in  Figure  1  in  Appendix  3.  The  conclusion  is  in  good  agreement  with 
experimental  results,  as  displayed  in  Fig.  6  of  Appendix  3.  The  characteristic  change  of  circular 
polarization  of  backscattering  light  with  size  of  scatterers  may  provide  a  new  method  to  directly 
determine  microscopic  size  of  scatterers  which  is  useful  for  highlighting  tumors  from  nonnal 
tissues. 


5.5.  Determination  of  light  absorption,  scattering  characteristics  and  anisotropy  factor 

We  are  developing  a  technique  to  determine  scattering  coefficient  (//v),  absorption  coefficient 
(jua),  and  anisotropy  factor  (g)  of  breast  tissue  specimens  using  backscattering  of  light  of 
different  wavelengths  (TO  2,  Task  4).  The  approach  involves  illumination  of  the  surface  of  the 
specimen  with  circularly  polarized  light  at  an  oblique  angle  of  incidence  (-70°)  and  measurement 
of  the  backscattered  co-polarized  and  cross-polarized  light.  From  the  backscattered  diffuse 
reflectance  data  the  absorption  coefficient  (//„),  and  the  reduced  scattering  coefficient  (ps )  are 
determined  at  different  wavelengths.  The  wavelength  dependence  of  the  reduced  scattering 
coefficient,  //v  is  used  to  extract  the  fractal  dimension  (Dj)  of  the  specimen.  The  depolarized 
diffuse  reflection  data  and  these  three  parameters  are  then  used  to  extract  the  anisotropy  factor 
(g).  The  g-factor  provides  a  useful  probe  of  the  size  of  the  scatterers  and  serves  as  a  fingerprint 
for  discriminating  between  normal  and  cancerous  tissues.  The  approach  may  prove  useful  in 
lumpectomy,  where  the  surgeon  looks  for  the  margins  of  the  tumor  with  high  accuracy. 

Details  of  the  experimental  method  and  approach  to  extract  key  parameters  are  presented  in 
Appendix  4\  “ Determination  of  light  absorption,  scattering  and  anisotropy  factor  of  a  highly 
scattering  medium  using  backscattered  circular  polarized  light,"  Proceedings  of  SPIE  Vol.  6435, 
64350J  (2007). 

6.  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Demonstrated  the  efficacy  of  Optical  Tomography  using  Independent  Component  Analysis 
(OPTICA)  approach  for  detection,  3-D  localization,  and  cross  section  imaging  of  a 
fluorescent  target  in  breast  tissue  specimens  and  model  medium.  (Appendix  1 ). 

•  The  location  of  a  tumor  inside  a  realistic  breast  model  composed  of  excised  breast  tissues 
was  determined  with  millimeter  accuracy  (Appendix  2). 

•  Developed  a  forward  model  for  detennining  the  backscattered  polarized  light  intensity 
distribution  (Appendix  3). 

•  Developing  a  technique  to  detennine  scattering  coefficient  (//,),  absorption  coefficient  (jua), 
and  anisotropy  factor  (g)  of  breast  tissue  specimens  that  has  potential  application  in  optical 
biopsy,  that  is,  for  detennining  if  the  tissue  specimen  is  nonnal  or  cancerous  (Appendix  4). 
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scattering  and  anisotropy  factor  of  a  highly  scattering  medium  using  backscattered  circular 
polarized  light.”  Paper  6935-18  presented  at  Conference  6435:  Optical  Interactions  with 
Tissue  and  Cells  XVIII  of  SPIE’s  BiOS  2007/  Photonics  West,  January  20-25,  2007,  San 
Jose,  California. 

3.  M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen  and  R.  R.  Alfano,  “Optical  high  resolution  cross  section 
imaging  of  human  breast  using  independent  component  analysis.”  Presented  at  the  5th 
Inter-Institute  Workshop  on  Optical  Diagnostic  Imaging  from  Bench  to  Bedside  at  the 
National  Institutes  of  Health,  25-27  September  2006,  Bethesda,  Maryland,  USA, 

8.  CONCLUSION 

The  work  carried  out  during  this  reporting  period  shows  the  potential  for  detection  and  three- 
dimensional  localization  of  targets  (including  tumor)  within  a  turbid  medium  (including  model 
breasts  fonned  with  ex  vivo  tissues)  with  significant  accuracy  based  on  the  differences  in  the 
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scattering,  absorption,  and  fluorescence  characteristics  of  the  target  ( e.g tumor)  and  the 

intervening  medium  (e.g.,  normal  breast  tissue). 

“So  What  Section” 

•  The  completion  of  the  Ph.  D  thesis  by  the  graduate  student  supported  by  the  grant  is 
important  since  it  adds  another  scientist  to  the  ongoing  struggle  against  cancer. 

•  A  recent  study  involving  35,3 19  patients  underscores  the  influence  of  primary  tumor  location 
on  breast  cancer  prognosis,  and  makes  it  imperative  that  breast  cancer  detection  modalities 
obtain  three  dimensional  (3-D)  location  of  the  tumor  relative  to  the  axilla.8  The  current  work 
is  an  important  development  in  obtaining  noninvasive  3-D  location  of  a  tumor  within  the 
breast. 

•  The  applicability  of  OPTICA  for  scattering,  absorptive,  and  fluorescent  targets  makes  it 
versatile  since  all  three  phenomena  may  be  used  for  contrast  enhancement  between  the  tumor 
and  normal  breast  tissues.  Fluorescence-based  detection  may  require  use  of  contrast  agents, 
such  as,  molecular  beacons. 

•  Three-dimensional  target  localization  will  enable  closer  probing  of  a  smaller  volume  around 
of  the  target  providing  more  details  since  smaller  pixel  size  could  be  used  without  increasing 
the  computation  time  (as  a  smaller  volume  will  be  probed). 

•  Further  development  of  the  OPTICA,  and  the  development  of  the  oblique  angle 
backscattering  have  potential  application  in  lumpectomy,  where  the  surgeon  looks  for  the 
margins  of  the  tumor  with  high  accuracy. 

•  What  is  more  important,  the  work  accomplished  so  far  prepares  us  for  the  next  step,  in  vivo 
optical  breast  imaging  involving  volunteers. 
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Localization  and  cross  section  reconstruction  of  fluorescent  targets  in  ex 
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An  information  theory  based  approach  is  introduced  for  the  detection  and  three-dimensional 
localization  of  fluorescent  targets  embedded  in  a  turbid  medium.  The  approach  uses  multisource 
illumination  of  the  medium,  multidetector  transillumination  acquisition  of  fluorescence  signal, 
independent  component  analysis  for  target  detection  and  localization,  and  a  back-projection 
algorithm  for  reconstruction  of  target  cross  section.  The  efficacy  of  the  approach  is  demonstrated  by 
locating  and  estimating  the  cross  section  of  a  fluorescent  target  embedded  in  a  100X100 
X  26  mm3  ex  vivo  human  breast  tissue  specimen  and  in  a  tissue-simulating  turbid  medium  of 
thickness  50  times  the  transport  mean  free  path.  ©  2006  American  Institute  of  Physics. 

[DOI:  10.1063/1.2356024] 


Development  of  optical  fluorescence  tomography  is  of 
interest  for  cancer  detection  and  diagnosis. 1-4  The  recent  ad¬ 
vent  of  fluorescence  beacons  and  contrast  agents  with  the 
ability  to  attach  themselves  to  desired  abnormal  cells  or 
organelles5  enhances  the  target  specificity  and  diagnostic  po¬ 
tential  of  fluorescence  based  optical  imaging  approaches. 
Determination  of  accurate  location  of  targets  embedded  in  a 
turbid  medium  remains  a  formidable  task  for  both  direct  im¬ 
aging  and  inverse  image  reconstruction  approaches  and  for 
all  types  of  targets — fluorescent,  absorptive,  and 
scattering. 6-11  We  have  introduced  an  approach,  known  as 
optical  tomography  using  independent  component  analysis 
(OPTICA),  that  enables  more  accurate  determination  of  tar¬ 
get  location  in  turbid  media  than  existing  optical 
methods.1-'13  It  uses  multisource  optical  probing  of  the 
sample,  multidetector  acquisition  of  transmitted  signal,  and 
independent  component  analysis14  of  information  theory  to 
analyze  the  signal  for  three-dimensional  (3D)  localization 
and  characterization  of  the  target.  The  efficacy  of  OPTICA 
for  locating  small  absorptive  and  scattering  targets  embedded 
in  turbid  media  has  been  demonstrated.11'15  In  this  letter, 
OPTICA  formalism  is  extended  for  fluorescent  targets,  and 
its  efficacy  is  tested  by  obtaining  3D  location  and  cross  sec¬ 
tional  image  of  a  fluorescent  target  in  an  ex  vivo  human 
breast  tissue  sample  and  in  a  tissue-simulating  turbid  me¬ 
dium  from  experimental  data. 

The  theoretical  formalism  for  fluorescence  OPTICA  is 
based  on  the  premise  that  the  spatial  distribution  of  the  light 
intensity  at  the  exit  boundary  of  the  medium  is  a  weighted 
mixture  of  signals  arriving  from  the  fluorescent  targets 
(“sources”)  embedded  in  the  medium: 

4rd,rs)  =  2  a/rjs/rj ,  (1) 

1 

where  afrf  is  the  mixing  vector  and  a/rd)s/rs)  represents 
the  contribution  of  the  /'th  target  to  the  fluorescent  signal  at 
the  detector  plane  for  illumination  of  the  source  plane  at  iy 
Independent  component  analysis  assumes  that  these  fluores¬ 
cent  sources  are  independent,  treats  the  problem  of  detecting 
those  as  a  source  separation  problem,  and  obtains  cy-  and  Sj 
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by  seeking  maximal  mutual  independence  between  s/s.12'14 

Light  propagation  in  a  highly  scattering  medium  with  an 
embedded  fluorescent  target  excited  by  an  external  light 
source  is  approximately  described  by  coupled  diffusion 
equations  at  the  excitation  and  emission  wavelengths.16'17 
The  fluorescence  signal  Um{ rd,rs,(o)  can  be  expressed  in 
terms  of  the  two  Green’s  functions  G/r.iyw)  and 
G,„(rrf,r.oj)  describing  light  propagation  from  the  excitation 
source  at  r5  to  the  fluorescent  target  at  r  at  the  excitation 
wavelength  kx  and  the  light  propagation  from  the  target  to 
the  detector  at  rd  at  the  emission  wavelength  X.,,,.  In  this  letter 
we  report  on  slab  samples.  Consequently,  Green’s  functions 
for  a  slab  geometry18  under  the  diffusion  approximation  of 
the  equation  of  radiative  transfer  are  used. 

Assuming  that  the  /th  fluorescent  target  is  contained  in 
volume  V j  centered  at  iy  the  fluorescence  signal  under  exci¬ 
tation  by  a  unit  point  source  at  rv  is  given  by 

1 

U„,( rd,rs,to)  =  2  Gm(rd,  iy  co)// w)  Gx{ iy  iy  w) ,  (2) 

j=i 

where  //to)  =  y(r/c  Vy-/[  1  -/coAiy)]  represents  the  fluores¬ 
cence  signal  strength  of  the  y'th  target,  y  is  the  fluorescent 
yield,  c  is  the  speed  of  light  in  medium,  t  is  the  fluorescence 
lifetime,  and  to  is  the  angular  modulation  frequency  of  the 
excitation  light  intensity. 

Both  the  location  and  the  strength  of  the  j  th  target  can  be 
computed  by  using  a  least  square  fitting  procedure: 

min  j  2  [</“ '/(g)  -  G/iyr,.)]2  +  2  rd) 

rj’aj’Pj  [  rs  rd 

-Gm(rd,r;)]2j,  (3) 

where  ,s/rs)  «  G/ry ,  rs. ,  of  and  a/iy)  ^  Gm(rd ,  r;- ,  at) .  The  fit¬ 
ting  yields  the  location  r;  of  and  the  two  scaling  constants  cy¬ 
an  d  [3j  for  the  /th  target.  The  fluorescence  strength  then  is 

fj=aA 

The  size  and  shape  of  the  /th  target  can  be  estimated 
from  a  back  projection  of  Um.(rd,rs,  w)  from  the  detection 
plane  onto  the  “target  plane,”  i.e.,  Z=Zj  plane.  The  fluores¬ 
cence  signal  from  the  /th  target  is  approximated  by 
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Um.{  rd 


I 


r  s,co)=  Gm(pd- p,(o)Xj{p)Gx(p-  ps,co)dp. 


Scanning 


(4) 


where  the  integration  is  over  the  z=Zj  plane  and  pd  and  ps  are 
the  lateral  coordinates  of  the  detector  and  the  source,  respec¬ 
tively.  In  the  Fourier  space,  X;(q)  can  be  obtained  from  Eq. 
(4)  as 


*i(q) : 


Gm((l  -  qJ5  w)G*(qj;to)  ’ 


(5) 


where  q  and  qs  are  the  spatial  frequencies  on  the  x-y  plane 
and  denotes  complex  conjugate.  The  inverse  Fourier  trans¬ 
forms  of  Xy(q)  yield  the  cross  sectional  image  of  the  fluores¬ 
cent  target  in  the  z=Zj  plane.  The  full  width  at  half  maximum 
(FWHM)  of  the  cross  sectional  image  provides  an  estimate 
of  the  real  target  size.  The  cross  sectional  image  is  a  two- 
dimensional  (2D)  distribution  of  the  fluorescent  object  at  the 
target  plane. 

The  experimental  arrangement  to  test  the  fluorescence 
OPTICA  approach  is  shown  schematically  in  Fig.  1(a).  The 
fluorescent  target  was  a  glass  sphere  (outer  diameter  of 
~4  mm  and  inner  diameter  of  3.2  mm)  filled  with  a  solution 
of  indocyanine  green  (ICG)  in  water.  The  water  solution  of 
ICG  had  an  absorption  coefficient  of  11.5  cm-1  at  784  nm 
and  fluoresced  over  the  790-966  nm  spectral  range  with  a 
peak  at  825  nm.  The  sphere  was  placed  inside  a  100X  100 
X  26  mm3  slab  of  ex  vivo  female  breast  tissue  specimen  con¬ 
tained  in  a  transparent  plastic  box.  One  of  the  sides  of  the 
box  could  be  moved  to  uniformly  compress  and  hold  the 
tissue  specimen  in  position.  The  breast  tissue  specimen  con¬ 
sisted  primarily  of  adipose  tissue  with  small  amount  of  skin 
and  was  provided  to  us  by  the  National  Disease  Research 
Interchange  under  an  Internal  Review  Board  approval  at  the 
City  College  of  New  York.  In  another  experiment  to  test  if 
the  approach  could  locate  a  target  inside  an  average-size  hu¬ 
man  female  breast,  the  sphere  was  placed  inside  a  250 
X  250  X  50  mm3  transparent  plastic  cell  filled  with  Intralipid- 
10%  suspension  in  water.  The  concentration  of  Intralipid- 
10%  was  adjusted  to  provide  transport  lengths  /,  of  —  1  mm 
at  784  nm  and  1.1  mm  at  830  nm,  which  are  within  the  re¬ 
ported  range  of  values  for  healthy  human  breast  tissues.19 

A  200  pm  optical  fiber  delivered  a  784  nm,  300  mW 
continuous-wave  beam  from  a  diode  laser  for  sample  illumi¬ 
nation.  The  beam  was  collimated  to  a  1  mm  spot  onto  the 
entrance  face  (henceforth  referred  to  as  the  “source  plane”) 
of  the  slab  samples.  Multiple  source  illumination  was  real¬ 
ized  in  practice  by  step  scanning  the  slab  samples  across  the 
laser  beam  in  an  x-y  array  of  grid  points  using  a  computer 
controlled  translation  stage.  The  x-y  array  was  18  X  15  for 
the  tissue  sample  and  10  X  10  for  the  Intralipid- 10%  sample, 
with  a  step  size  of  2.5  mm  for  both  the  cases.  The  signal 
from  the  opposite  face  of  the  sample  (henceforth  referred  to 
as  the  “detection  plane”)  was  passed  through  a  narrow-band 
interference  filter  centered  at  830  nm  (FWHM  —10  nm,  50% 
transmission)  to  block  the  scattered  784  nm  pump  light.  A 
camera  lens  collected  the  fraction  of  filtered  fluorescence 
around  830  nm  and  projected  it  onto  the  sensing  element  of  a 
cooled  16  bit,  1024  X  1024  pixel  charged  couple  device 
(CCD)  camera.  Each  illuminated  pixel  of  the  CCD  camera 
could  be  regarded  as  a  detector.  For  illumination  of  every 

scanned  point  on  the  source  plane,  the  CCD  camera  recorded 
.  .  -  - -t0  fr-  . . 


l-ll-El 


Cropped  Burned 

image  image 


FIG.  1.  (a)  Schematic  diagram  of  the  experimental  arrangement  (CCD 
= charge  coupled  device  camera,  F= filter,  PC = personal  computer),  (b)  An 
expanded  view  of  the  sample  input  (source)  plane  and  exit  (detection)  plane 
showing  the  grid  points  in  the  x-y  plane,  (c)  A  photograph  of  the  exit  face  of 
the  tissue  sample,  (d)  Absorption  and  fluorescence  spectra  of  the  ICG  dye 
solution  in  water,  (e)  A  typical  raw  CCD  image  of  the  detection  plane  and 
how  it  is  cropped  and  binned  for  analysis. 


the  diffusely  transmitted  2D  fluorescence  intensity  pattern  on 
the  detection  plane.  A  second  scan  to  estimate  the  average 
value  of  k=  y3fia/LLrs  (where  fia  and  fis  are  the  absorption  and 
scattering  coefficients,  respectively)  of  the  background  me¬ 
dium  used  a  short-pass  filter  at  750  nm  to  block  the  fluores- 


30 

x?y  (mm) 

FIG.  2.  OPTICA  generated  fluorescence  intensity  distribution  on  (a)  the 
detector  and  (b)  the  source  planes,  [(c)  and  (d)]  Green’s  function  fits  to  the 
horizontal  spatial  profiles  through  the  centers  of  intensity  distributions  in  (a) 
and  (b),  respectively,  (e)  Cross  sectional  image  at  the  z=Zj  plane,  (f)  Spatial 
profiles  of  the  cross  sectional  image  along  the  x  and  y  directions  shown  by 
the  white  lines  in  (e). 
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TABLE  I.  Comparison  of  known  positions  and  sizes  of  the  target  with  OPTICA-determined  values. 


Known  position 

Observed  position 

Profile  FWHMa 

Samples 

(x,  y,  z)  (mm) 

(x,  y,  z)  (mm) 

(Ax,  Ay)  (mm) 

Breast  tissue  ex  vivo 

(30,  27,  15) 

(30,  27,  16) 

(4,  4) 

Intralipid- 10%  in  water 

(115,  115,  15) 

(115,  115,  16) 

(12,  12) 

“The  FWHM  values  should  be  compared  with  the  inner  diameter  of  3.2  mm  of  the  sphere. 


cence  signal  and  admit  a  very  small  fraction  of  the  transmit¬ 
ted  excitation  light. 

Figure  1(b)  shows  a  photograph  of  the  exit  surface  of  the 
tissue  sample.  Figure  1(c)  presents  a  typical  2D  raw  image 
(fluorescence  intensity  distribution)  of  the  detection  plane 
recorded  by  the  CCD  camera  for  illumination  of  a  grid  point 
in  the  source  plane.  Each  raw  image  is  then  cropped  to  select 
out  the  information-rich  region  and  binned  to  enhance  the 
signal-to-noise  ratio.  All  the  binned  images  corresponding  to 
illumination  of  the  grid  points  in  sequence  are  then  stacked 
and  used  as  input  for  independent  component  analysis. 

Figures  2(a)  and  2(b)  present  the  OPTICA  generated  in¬ 
tensity  distributions  on  the  detector  and  source  planes,  re¬ 
spectively,  of  the  tissue  sample.  The  corresponding  Green’s 
function  fits  to  horizontal  profile  through  the  center  of  inten¬ 
sity  distributions  are  shown  in  Figs.  2(c)  and  2(d),  respec¬ 
tively.  The  (x,y  ,z)  location  of  the  target  center  obtained  from 
this  analysis  is  displayed  in  Table  I.  The  lateral  positions 
( x,y )  agree  completely,  and  the  axial  position  agrees  within 
1  mm  of  the  known  values.  The  cross  sectional  image  ob¬ 
tained  using  back-projection  Fourier  transforms  in  Eq.  (5),  is 
shown  in  Fig.  2(e).  Figure  2(f)  shows  intensity  profiles  of  the 
cross  sectional  image  along  the  x  and  y  directions  shown  by 
the  white  dotted  lines  in  Fig.  2(e).  The  FWHM  values  of 
intensity  profiles  of  the  cross  sectional  image,  presented  in 
column  4  of  Table  I,  are  4  mm,  in  good  agreement  with  the 
known  3.2  mm  inner  diameter  of  the  fluorescent  sphere. 

The  results  of  measurements  on  the  same  fluorescent 
sphere  embedded  in  Intralipid- 10%  suspension  in  water  are 
also  summarized  in  Table  I.  Again  the  lateral  positions  are  in 
exact  agreement,  while  the  axial  position  agrees  within 
1  mm.  However,  the  FWHM  values  of  the  intensity  profiles 
of  the  cross  sectional  image  are  estimated  to  be  12  mm  each, 
which  is  3.7  times  the  3.2  mm  inner  diameter  of  the  sphere. 
We  attribute  this  difference  primarily  to  the  much  larger 
thickness  (50/,)  of  the  Intralipid- 10%  sample  than  that  of  the 
tissue  sample  (26/,).  Spatial  resolution  decreases  rapidly  with 
sample  thickness  (for  the  same  axial  target  position)  because 
of  the  diffuse  nature  of  light  propagation  through  a  highly 
scattering  medium.  Since  the  transport  lengths  (/,)  were  com¬ 
parable  and  the  target  was  located  15  mm  from  the  source 
plane  for  both  the  samples,  the  target-detector  plane  distance 
in  Intralipid- 10%  sample  was  more  than  three  times  larger 
than  that  in  the  tissue  sample,  which  in  turn  led  to  a  much 
higher  estimate  for  the  target  size.  As  the  incident  beam 
power  was  limited  to  300  mW  for  both  the  samples,  higher 
optical  thickness  resulted  in  much  lower  fluorescence  signal- 
to-noise  ratio  for  the  Intralipid- 10%  sample,  which  in  turn 
contributed  to  a  higher  estimate  of  the  target  cross  section. 
Substantial  improvements  in  target  size  estimate  are  ex¬ 
pected  if  an  intensified  CCD  camera  for  higher  sensitivity, 
and  more  incident  beam  powers  are  used. 


It  should  be  noted  that  multisource  measurements, 
though  crucial  for  detection  of  multiple  targets,  might  not  be 
needed  for  a  single  target,  since  the  position  of  the  target 
determines  the  normalized  fluorescence  intensity  distribu¬ 
tion.  However,  multisource  measurements  help  reduce  the 
deleterious  effects  of  noise  and  inaccuracy  in  the  estimation 
of  optical  properties  of  the  medium  and  provide  better  as¬ 
sessment  of  3D  location  even  in  the  case  of  a  single  target. 

The  simpler  but  important  case  involving  a  single  target 
(e.g.,  a  tumor  in  the  breast)  was  chosen  for  experimental 
demonstration.  In  simulation  the  approach  could  locate  three 
targets  with  one  target  —20/,  behind  another  in  a  —50/,  thick 
turbid  medium.  Localization  of  even  more  targets  is  possible 
and  will  be  presented  along  with  experimental  results  in  a 
future  publication. 

In  summary,  OPTICA  can  provide  rather  accurate  3D 
location  of  a  fluorescent  target  in  a  turbid  medium  and  an 
estimate  of  its  cross  section. 

This  work  is  supported  in  part  by  PSC-CUNY,  U.S. 
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ABSTRACT 

Optical  imaging  using  independent  component  analysis  (OPTICA)  is  enhanced  to  provide  a  high  resolution  cross  section 
imaging  of  objects  in  a  turbid  medium  by  a  backprojection  technique.  The  performance  is  demonstrated  by  imaging  a 
human  breast  model  made  of  ex  vivo  human  breast  tissues.  Cancerous  site  of  5mm  size  is  detected  at  the  midplane  of  the 
33mm  thick  breast  model.  The  reconstructed  cross  section  image  compares  favorably  with  pathology  findings. 

Keywords:  optical  imaging,  independent  component  analysis,  cross  section,  human  breast 


1.  INTRODUCTION 

Biomedical  diffuse  optical  tomography  (DOT)  has  received  significant  attention  and  development  in  the  past  two  decades 
and  has  been  applied  to  image,  such  as,  human  breast,  infant  brain  and  joints.  The  state  of  art  diffuse  optical  tomography 
currently  relies  on  an  iterative  image  reconstruction  approach  to  construct  the  optical  properties  of  the  interior  of  the 
medium.  This  approach  is  computationally  demanding  yet  provides  limited  spatial  resolution. 1  Reconstruction  of  images 
with  adequate  spatial  resolution  and  accurate  localization  and  characterization  of  the  inhomogeneities  remain  a  formidable 
task. 

An  alternative  optical  imaging  approach  is  optical  imaging  using  independent  component  analysis  (OPTICA).  The 
principle  of  OPTICA  is  to  first  sort  out  the  signal  due  to  individual  inhomogeneity  and  then  perform  three  dimensional  (3D) 
localization  and  characterization  for  each  inhomogeneity,  using  a  massive  data  set  generated  by  a  multi-source  illumination 
and  multi -detector  signal  acquisition  scheme.  Independent  component  analysis  has  been  extensively  applied  in  fields  such 
as  Electroencephalogram  (EEG)  and  functional  MRI  (fMRI)  and  shown  effective  in  separating  signals  from  different  brain 
activity  centers.  Our  previous  work  has  shown  that  such  an  approach  significantly  improves  the  sensitivity  to  small/weak 
absorptive,  scattering  and/or  fluorescence  inhomogeneities  and  can  achieve  a  3D  localization  of  the  target  with  error  less 
than  1mm. 

In  this  paper,  we  report  on  a  back  projection  method  to  obtain  high-resolution  cross-section  images  of  targets  inside 
highly  scattering  media,  enhancing  OPTICA.  We  demonstrate  the  performance  of  OPTICA  by  imaging  a  42  x  30  x  33mm  3 
model  human  breast  assembled  using  ex  vivo  breast  tissues.  A  5  x  5  x  3mm3  tumor  was  placed  at  the  midplane.  OPTICA 
located  the  tumor  within  ~  1mm  of  its  known  position.  It  also  identified  a  10mm  low  scattering  structure  at  midplane,  and 
a  3mm  low  scattering  structure  near  the  surface.  Subsequent  pathological  analysis  confirmed  the  tumor,  and  identified  the 
other  two  structures  as  glandular  breast  tissues.  The  reconstructed  cross  section  images  and  optical  properties  of  the  targets 
were  consistent  with  pathology  findings. 
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2.  THEORETICAL  FORMALISM  OF  OPTICAL  IMAGING  USING  INDEPENDENT  COMPONENT 

ANALYSIS 

The  presence  of  optical  inhomogeneities  inside  a  turbid  medium  perturbs  the  intensity  of  emergent  light  from  the  medium 
under  illumination  by  a  probing  beam.  In  the  diffusion  approximation,  the  change  of  the  light  intensity  on  the  boundary  of 
the  specimen  due  to  absorptive  and  scattering  objects  (inhomogeneities)  can  be  written  as:  2,3 

-A/(rd,rs)  =  Jd*r5pa(r)cG(rd,r)G(r,rs)  +  J  d3r5D(r)cVrG(rd,  r)  •  VrG(r,  rs)  (1) 


to  the  first  order  of  Born  approximation4  when  illuminated  by  a  point  source  of  unit  power,  where  r s  and  r,/  are  the 
positions  of  the  probing  source  and  the  light  detector  on  the  boundary,  8p  a(r)  =  /ra(r)  —  pao  and  6D( r)  =  D( r)  —  Dq 
are  the  differences  in  absorption  coefficient  and  diffusion  coefficient,  respectively,  between  the  inhomogeneity  at  r  and  the 
background,  c  is  the  speed  of  light  in  the  medium,  and  G(r,  r ')  is  the  Green’s  function  describing  light  propagation  from 
r'  to  r  inside  the  background  turbid  medium  of  absorption  and  diffusion  coefficients  //  ao  and  Do- 

OPTICA  assumes  each  inhomogeneity  within  the  scattering  medium  to  be  a  virtual  source  and  express  the  change  of 
the  light  intensity  on  the  boundary  of  the  specimen  as: 

-AI(rd,rs)  =  '^2aj(rd)sj(rs)  (2) 

j 


where  Sj  (rs)  represents  the  j-th  inhomogeneity  illuminated  by  the  incident  wave  at  r  s,  and  aj  ( rd )  is  the  weighting  matrix 
describing  the  propagation  of  light  from  the  j-th  inhomogeneity  to  the  detector  at  r  d.  Each  absorptive  inhomogeneity 
contributes  one  term  in  Eq.  (2)  and  each  scattering  inhomogeneity  contributes  three  terms  in  Eq.  (2).  The  detected  change 
of  the  light  intensity  (2)  is  hence  a  linear  mixture  of  signals  where  a  j  and  Sj  can  now  be  interpreted  as  the  j-th  weighting 
matrix  and  virtual  source,  respectively.  Owing  to  the  statistical  independence  between  these  virtual  sources,  independent 
component  analysis  of  —A  I  will  yield  a  list  of  independent  components  and  recover  both  a  j  and  Sj.  The  number  of  the 
leading  independent  components  gives  the  number  of  objects .  The  location  and  characterization  of  the  j-th  target  is  obtained 
from  the  analysis  of  the  retrieved  independent  component  (s  j  and  a:l )  which  relates  directly  to  the  source-to-object  and 
object-to-detector  Green’s  functions  G(r_, ,  rs)  and  G(rd,  r j)  and  the  optical  property  of  the  target  where  r_,  is  the  position 
of  the  j-th  object. 5~8 

For  a  slab  geometry  investigated  here,  there  are  three  virtual  sources  of  specific  patterns  (one  centrosymmetric  and  two 
dumbbell-shaped)  associated  with  one  scattering  inhomogeneity,  whereas  only  one  centrosymmetric  virtual  source  is  asso¬ 
ciated  with  one  absorptive  inhomogeneity.  Among  the  three  virtual  sources  associated  with  one  scattering  inhomogeneity, 
the  centrosymmetric  virtual  source  is  the  easiest  to  detect.  The  centrosymmetric  virtual  source  and  the  corresponding 
weighting  matrix  are  Sj  oc  G(jj,rs)  and  aj  <x  G(rd,Vj),  and  Sj  ex  ^(rj,rs)  and  aj  ex  ^(r^re,),  respectively,  for 
absorptive  and  scattering  inhomogeneities.  A  simple  least  square  fitting  of  the  centrosymmetric  component,  such  as 


min 

r  j  ,ctj  ,/3j 


SOs)  -  G( r 


3  ’ 


)] 2  +  [Pj  1aj(rd)  -  G(rd,  ij)] ' 


(3) 


for  the  absorptive  object,  can  be  used  to  yield  the  3D  location  r  ,  and  the  strength  a:l  ;3:l  of  the  target.  When  a  prior 
knowledge  about  the  property  of  the  target  is  not  available,  Eq.  (3)  can  still  be  used  to  estimate  the  3D  location  of  the  target 
regardless  absorption  or  scattering  property  of  the  target  since  rs)  ~  — «G(rj,  rs)  and  ^f(rdl  rj)  ~  —kG(  rdl  r_,j 

where  k  =  sj (/zao  —  iu>/c)/ Do  chosen  to  have  a  nonnegative  real  part  with  ui  the  modulation  frequency  of  the  incident 
wave. 

The  signal  of  the  j-th  target  is  simply  given  by  —Alj  =  aj(rd)sj(rs).  On  the  other  hand,  the  centrosymmetric  signal 
of  the  j-th  target  can  be  approximately  expressed  as  a  double  convolution  in  the  form  of 

-AIj(rd,rs)=  f  G(pd- p,zd,Zj)Xj(p)G(p  -  ps,Zj,zs)dp  (4) 
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Figure  2.  Independent  light  intensity  distributions  on  the  detector  plane  obtained  by  OPTICA  for  the  cancer  site.  The  cancer  site  (C) 
was  found  at  a  depth  of  14.8mm. 


4.  RESULTS 

A  typical  raw  image  measured  on  the  detector  plane  at  one  scanning  position  is  shown  in  the  inset  of  Fig.  1.  The  average 
of  all  the  22  x  16  images  was  used  to  obtain  the  optical  property  of  the  background  medium.  The  absorption  coefficient 
of  the  background  medium  is  found  to  be  /ra  =  0. 0039mm-1  assuming  the  transport  mean  free  path  lt  =  1mm.  Each 
raw  image  is  first  cropped  to  retain  the  region  within  the  window  of  42  x  30mm  2  over  which  image  reconstruction  would 
be  performed,  and  binned  to  enhance  the  signal-to-noise  ratio.  All  the  binned  images  corresponding  to  illumination  of  the 
grid  points  in  sequence  were  then  stacked,  and  used  as  input  for  independent  component  analysis. 

The  independent  light  intensity  distributions  obtained  by  OPTICA  is  displayed  in  Fig.  2.  The  cancer  site  was  found  at 
the  depth  of  14.8mm  from  fitting  the  independent  light  intensity  distribution  to  the  background  Green’s  function. 

The  back-projection  of  the  detected  signal  for  the  cancer  site  (C)  produced  its  cross  section  images  at  the  depth  of 
14.8mm  in  Fig.  3.  The  corresponding  profile  along  the  horizontal  direction  is  displayed  on  the  right  for  the  target.  The 
cancer  site  is  estimated  to  be  11mm  in  size  and  it  scatters  light  more  strongly.  Representative  histology  of  the  site  is 
displayed  as  the  inset. 

The  investigated  ex  vivo  breast  sample  contained  minimal  amount  of  blood  and  hence  the  reconstructed  images  is  for 
the  scattering  property  of  the  sample.  The  change  of  /j for  the  targets  can  further  be  estimated  from  the  reconstructed 
independent  component  (sy  and  a:i )  for  the  site.  It  was  found  site  C  has  enhanced  scattering  compared  to  the  background 
(mainly  adipose  tissue).  Subsequent  pathological  analysis  confirmed  the  site  C  as  tumor. 

5.  DISCUSSION 

The  nature  of  the  inhomogeneity  (either  absorptive  or  scattering  or  mixed)  can  be  discovered  by  OPTICA  with  CW  mea¬ 
surement  when  the  signal  to  noise  ratio  is  high. 6’ '  When  the  signal  to  noise  ratio  is  not  favorable,  the  recovered  independent 
component  will  be  due  to  both  absorption  and  scattering  perturbations  at  the  site  of  the  inhomogeneity.  The  strength  of  the 
inhomogeneity  will  be  proportional  to  8fia  +  K25D  =  8/j,a  +  (/ iao  —iuj/c)8D/ Dq-  This  provides  a  very  simple  approach  to 
discriminate  absorption  vs  scattering  if  multiple  modulation  frequencies  ui  are  available.  The  capability  of  OPTICA  sepa¬ 
rating  absorption  from  scattering  inhomogeneities  can  be  significantly  improved  with  a  time-domain  or  frequency  domain 
measurement. 
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Analytical  cumulant  solution  of  the  vector  radiative  transfer  equation  investigates  backscattering 
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The  backscattering  of  circularly  polarized  light  pulses  from  an  infinite  uniform  scattering  medium  is  studied 
as  a  function  of  helicity  of  the  incident  light  and  size  of  scatterers  in  the  medium.  The  approach  considers  a 
polarized  short  pulse  of  light  incident  on  the  scattering  medium,  and  uses  an  analytical  cumulant  solution  of  the 
vector  radiative  transfer  equation  with  the  phase  matrix  obtained  from  the  Mie  theory  to  calculate  the  temporal 
profile  of  scattered  polarized  photons  for  any  position  and  any  angle  of  detection.  The  general  expression  for 
the  scattered  photon  distribution  function  is  an  expansion  in  spatial  cumulants  up  to  an  arbitrary  high  order. 

Truncating  the  expansion  at  the  second-order  cumulant,  a  Gaussian  analytical  approximate  expression  for  the 
temporal  profile  of  scattered  polarized  photons  is  obtained,  whose  average  center  position  and  half  width  are 
always  exact.  The  components  of  scattered  light  copolarized  and  cross  polarized  with  that  of  the  incident  light 
can  be  calculated  and  used  for  determining  the  degree  of  polarization  of  the  scattered  light.  The  results  show 
that  circularly  polarized  light  of  the  same  helicity  dominates  the  backscattered  signal  when  scatterer  size  is 
larger  than  the  wavelength  of  light.  For  the  scatterers  smaller  than  the  wavelength,  the  light  of  opposite  helicity 
makes  the  dominant  contribution  to  the  backscattered  signal.  The  theoretical  estimates  are  in  good  agreement 
with  our  experimental  results. 
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I.  INTRODUCTION 

The  recent  increased  interest  in  the  study  [1-12]  of  polar¬ 
ized  light  propagation  through  turbid  media  derives  from  a 
variety  of  practical  applications.  Polarized  light  multiple 
scattered  from  a  turbid  medium  carries  information  about  the 
medium  interior  [13]  such  as,  particle  size  and  concentration, 
size  distribution,  and  refractive  index  variation,  which  can  be 
useful  for  biomedical  imaging  [13-15],  endoscopic  evalua¬ 
tion  of  biological  tissues  [16],  investigation  of  biological  cell 
differentiation  [17,18],  flow  cytometry  [18,19],  lidar-based 
remote  sensing  of  the  atmospheric  cloud,  aerosol  fog  and 
smog  [20,21],  and  imaging  of  targets  in  shallow  coastal  wa¬ 
ter  [22].  Polarization-sensitive  imaging  has  been  demon¬ 
strated  to  enhance  image  contrast  of  targets  embedded  in 
turbid  media  [14,23]. 

In  particular,  it  has  been  observed  that  for  circularly  po¬ 
larized  light,  randomization  of  its  polarization  requires  more 
scattering  events  than  the  randomization  of  its  direction 
[23-27,6,9].  An  important  consequence  of  this  “polarization 
memory”  [25]  is  higher  resolution  and  contrast  of  transillu¬ 
mination  images  recorded  with  copolarized  light  than  those 
recorded  with  unpolarized  or  cross-polarized  light.  Another 
interesting  feature  is  the  dependence  of  the  intensity  of  the 
backscattered  copolarized  and  cross-polarized  light  on  the 
size  of  the  scatterers  in  the  medium  [23].  The  backscattered 
light  is  dominated  by  light  of  helicity  opposite  to  that  of  the 
incident  light  if  the  size  of  the  scatterers  is  smaller  than  the 
wavelength  of  light.  The  light  of  the  same  helicity  dominates 
when  scatterers  are  larger  than  the  light  wavelength.  These 
results  have  important  consequences  in  the  design  consider¬ 
ations  of  imaging  systems  using  backscattered  light  such  as, 
endoscopic  systems.  All  these  developments  make  it  impera¬ 


tive  that  theoretical  formalisms  be  developed  to  provide  a 
quantitative  explanation  of  the  experimental  results,  and 
guidance  for  the  design  of  imaging  systems  using  polarized 
light. 

Theoretical  studies  of  polarized  light  propagation  through 
a  scattering  medium  commonly  starts  with  the  vector  radia¬ 
tive  transfer  equation  (VRTE)  [1],  and  because  of  the  inher¬ 
ent  complexity  of  the  problem,  a  limited  number  of  studies 
are  available  to  date.  A  majority  of  these  studies  focused  on 
developing  numerical  solutions  of  the  VRTE.  Ishimaru  and 
co-workers  first  investigated  the  propagation  of  continuous 
wave  (CW),  polarized  wave  normally  incident  on  a  medium 
with  discrete  scatterers  using  discrete  ordinate  methods  to 
solve  VRTE  [2,3].  More  recently,  they  extended  the  formal¬ 
ism  to  numerically  solve  the  time-dependent  VRTE  to  ac¬ 
count  for  the  propagation  of  polarized  pulses  [4].  Kim  and 
Moscoso  studied  the  propagation  of  a  linearly  polarized  or 
circularly  polarized  continuous  plane  wave  normally  incident 
on  a  plane-parallel  medium  containing  a  random  distribution 
of  dielectric  spheres,  and  numerically  solved  the  one¬ 
dimensional  VRTE  using  a  Chebyshev  spectral  method  [5]. 
The  approach  was  subsequently  modified  to  study  the  back- 
scattering  of  circularly  polarized  plane-wave  pulses  [6].  Val- 
lion  et  al.  [7]  used  a  vector  Monte  Carlo  method  to  study 
polarized  light  transport  through  a  semitransparent  medium 
filled  with  scattering  particles.  Jiang  et  al.  [8]  developed  a 
multilayer  radiative  transfer  model  that  takes  into  account 
the  state  of  polarization  using  the  doubling-adding  method 
for  passive  atmospheric  remote  sensing  applications.  More 
recently,  Sakami  and  Dogariu  [9]  investigated  the  propaga¬ 
tion  of  a  polarized  plane-wave  pulse  incident  at  any  angle  in 
a  random  medium  using  the  discrete-ordinates  method  to 
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solve  the  VRTE,  and  used  an  angular  analysis  to  obtain  de¬ 
tails  of  the  polarization  flip  of  circular  polarized  light.  The 
problem  of  light  scattering  and  radiative  transfer  in  scattering 
media  and  its  application  to  geophysical  optics,  image  trans¬ 
fer,  remote  sensing,  and  inverse  problems  are  discussed  in 
the  monograph  by  Kokhanovsky  [10]. 

In  this  paper,  we  use  the  analytical  solution  of  the  time- 
dependent  VRTE  to  compute  backscattering  of  circularly  po¬ 
larized  light  in  an  infinite  uniform  scattering  medium.  It 
builds  on  the  analytical  cumulant  approach  for  solving  the 
time-dependent  radiative  transfer  equation  in  an  infinite  uni¬ 
form  medium  that  we  developed  [28-30].  In  this  approach  an 
arbitrary  phase  function,  as  long  as  it  is  a  function  of  scat¬ 
tering  angle,  can  be  treated.  We  have  derived  the  exact  ex¬ 
pression  of  spatial  cumulants  of  photon  distribution  at  any 
angle  and  any  time,  up  to  an  arbitrary  high  order.  The  exact 
first  spatial  cumulant  represents  the  center  of  the  distribution 
and  the  exact  second  spatial  cumulant  represents  the  width  of 
the  distribution.  The  method  of  the  analytical  cumulant  solu¬ 
tion  for  the  scalar  radiative  transfer  equation  can  be  extended 
to  the  vector  (polarized)  case  [31].  In  the  analytical  solution 
of  VRTE,  a  transform  to  the  circular  representation  from  the 
Stokes  representation  of  polarization  is  introduced,  and  the 
generalized  spherical  functions  are  used  for  angular  expan¬ 
sion,  instead  of  the  standard  spherical  functions  for  solving 
the  scalar  radiative  transfer  equation.  The  expressions  for  the 
exact  first  and  second  spatial  cumulants  of  polarized  compo¬ 
nents  as  functions  of  angle  s  and  time  t  is  derived.  Then, 
intensity  of  the  polarized  light  is  presented  by  a  Gaussian 
spatial  distribution.  We  apply  this  method  to  calculate  the 
backscattering  of  circular  polarized  light  in  a  scattering  me¬ 
dium.  The  results  show  that  circularly  polarized  light  of  the 
same  helicity  dominates  the  backscattered  signal  when  scat- 
terer  size  is  larger  than  the  wavelength  of  light.  For  the  scat¬ 
tered  smaller  than  the  wavelength,  light  of  opposite  helicity 
makes  the  dominant  contribution  to  the  backscattered  signal. 
The  theoretical  conclusion  is  in  good  agreement  with  our 
experimental  results  [23]. 

The  paper  is  organized  as  follows.  In  Sec.  II  we  review 
the  analytical  solution  of  VRTE.  In  Sec.  Ill,  the  numerical 
results  are  presented  and  are  compared  with  experiments. 
Section  IV  is  devoted  to  discussion  and  conclusion. 


(9I(r,s,f)/*  + cs  •  VrI(r,s,/)  +  /z„I(r,s,t) 

=  /jls  f  P(s,s')[I(r,s',t)  -  I(r,s,r)]afs' 

+ 1(0)  <5(r  -  r0)  <5(s  -  s0)  <5(t  -  0) ,  (1) 

where  vector  I(r ,  s ,  t)  has  four  polarization  components,  c  is 
the  light  speed  in  the  medium,  fis  is  the  scattering  rate,  j±a  is 
the  absorption  rate,  and  P(s,s')  is  the  4X4  phase  matrix.  A 
meridian  plane  parallel  to  the  z  axis  and  the  light  direction  s, 
is  used  as  a  plane  of  reference  for  the  description  of  the 
polarization  state.  In  SP  the  components  Q ,  U  vary  with  the 
rotation  of  the  reference  plane  around  the  light  propagation 
direction.  With  a  rotation  of  reference  plane  through  an  angle 
a  3*  0  (in  the  counterclockwise  direction,  when  looking  in  the 
direction  of  propagation)  Isp  varies  as  (l')sr=lJa)lsl\ 
where  L (a)  is  given  by 

10  0  0 

0  cos  2a  sin2a  0 

L(or)  =  .  (2) 

0  -  sin  2  a  cos  2  a  0 

_0  0  0  1 

When  light  propagates  along  the  z  direction,  Q ,  U  vary  with 
the  change  of  the  azimuthal  angle  </>. 

The  phase  matrix  in  the  fixed  coordinates  is  given  by: 

P(s,s')=L(i7-x)P(cos0)L(-x'),  (3) 

where  ©  is  the  angle  between  light  rays  before  and  after 
scattering  in  the  scattering  plane,  and  the  matrices  L(-^') 
and  are  those  required  to  rotate  meridian  planes 

before  and  after  scattering  onto  a  local  scattering  plane.  The 
intrinsic  property  of  a  scattering  mechanism  is  described  by 
the  4X4  scattering  function  P(cos  ©),  which  is  assumed  to 
depend  only  on  cos  @  =  s  s'.  It  is  convenient  to  use  a  repre¬ 
sentation  of  the  polarized  light  in  which  L {a)  is  diagonal, 
rather  than  Eq.  (2).  A  circular  polarization  representation 
(CP),  first  proposed  by  Kuscer  and  Ribaric  [34]  and  pre¬ 
sented  by  Hovenier  and  van  der  Mee  [35],  is  ICP 
=  [/2,/0,/_0,/_2],  where  70=(7+V)/2,  /_„  =  (/- V)/2,  I2  =  (Q 
+  iU)/2,  and  I_2=(Q-iU)  12,  or  IC7  =TE7’,  with 


II.  THEORETICAL  FORMALISM 
A.  Polarized  components  of  light 

In  the  Stokes  polarization  representation  [32]  (SP),  the 
polarized  light  is  described  by  I?p=[7,  Q,  U,  V\.  The  polar¬ 
ized  components  of  distribution  we  will  study  are  the  parallel 
component  Iu  =  (I+Q)/2,  the  perpendicular  component  / 

=  (I-Q)/2,  the  right-circular  component  IR=(I+V)/2,  and 
the  left-circular  component  IL=(I-V)/2,  since  these  quanti¬ 
ties  are  positive  definite,  and  are  experimentally  measurable. 

The  vector  radiative  transfer  equation  for  the  polarized 
photon  distribution  function  I(r,s,r)  in  an  infinite  uniform 
medium,  illuminated  by  a  point  light  source  providing  short 
pulse,  I(0)<S(r-r0)<5(s-s0)$7-0)  is  given  by  [33] 


0 

1 

i 

0 

0 

1 

1 

0 

1 

1 

0 

0 

1 

and  T-1  = 

1 

0 

0 

1 

T  =  - 

2 

1 

0 

0 

- 1 

-  i 

0 

0 

i 

0 

1 

-  i 

0 

0 

1 

-  1 

0 

(4) 

In  CP  a  rotation  of  the  reference  plane  through  an  angle  a 
around  the  light  direction  causes  lftlP  to  be  multiplied  by 

exp(-(/«a'),  m  =  2,0,0,-2.  Notice  that  70  and  / . „  actually 

have  the  same  rotational  property.  For  the  phase  matrix, 
transform  between  two  representations  is  Pr/=TPV,T  . 

In  CP  it  is  convenient  to  expand  the  phase  matrix  Pcp 
using  the  generalized  spherical  functions  (GSF).  The  gener¬ 
alized  spherical  functions  P'mn(/3 )  are  related  to  the  d'nm(/3) 
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in  the  angular  momentum  theory  [36]  by  Plm  n(f3) 
=  (i)m~ndImn(/3) .  In  the  Appendix,  the  expressions  for  d!mn(/3) 
are  presented. 

B.  Phase  matrix  using  Mie  theory 

Elements  of  the  CP  phase  matrix  in  the  scattering  plane, 
P,pn{ cos  ©)  can  be  expanded  by  GSF; 

pmn( cos  ©)  =  p'mnplm,n(C0S  ©)  .  (5) 

4tt  i 


where  s  =  (#,</>),  and  m,«0=2,0,-0,-2,  Z^max(|m|,|«0|). 
The  second  term  in  Eq.  (8)  represents  the  ballistic  (unscat¬ 
tered)  component  of  photons.  Flmn(t)  in  Eq.  (8)  is  the  solu¬ 
tion  of  the  equation 

^Lo(0/*=En'„X„o(0,  (9) 

n 

under  the  initial  condition 

plmn0^  =  0)  =  Sm„o(2l  +  1)/4tt,  (10) 


with  m , /z  =  2 , 0 ,  —0 ,  —2,  and  Z>max(|m|,|n|).  The  coeffi¬ 
cients  p'mn  provide  a  description  of  the  scattering  mechanism. 
In  most  useful  cases,  the  coefficients  pmn  have  properties: 

(l)  p'mm  =  Pm-m  are  real;  (“)  P‘mn  =  Pnm  =  P-m-n’  ("0  P20 
=  [pn_()J  ( *  means  complex  conjugate).  Therefore,  for  each 
Z 3=  2,  there  are  six  independent  real  elements:  p[)0,  p22,  Po_0, 
p(_2,  Re[p!,0],  and  Im[p(,0].  For  Z=0, 1,  only  Pq0  and  Pq_0  are 
nonzero.  These  numerical  coefficients  can  be  calculated  us¬ 
ing  the  Mie  theory  for  a  spherical  particle.  We  have  p'00 
=  (a[  +  a!4)/2,  p'0_0=(a[-a'4)/2,  p22  =  (a2  +  a3)/2,  p2_2  =  (a2 
-a'3)/2,  Re[p20]  = /3j ,  Im[p(,0]  = /32,  where  the  coefficients 
a\  2  3  4  and  f3\  2  are  formulated  in  Eqs.  (4.81)-(4.86)  in  Ref. 

[37]-' 


C.  Photon  angular  distribution  F(s,s0,f) 

When  we  make  a  spherical  harmonics  expansion  of  Eq. 
(1),  the  difficulty  is  that  the  term  cs- VrI(r,s,t)  couples 
equations  of  different  harmonics.  We  first  study  the  photon 
distribution  in  the  angular  space  in  CP  initially  incident  along 

s0, 


F(s,s0,f)  =  J  I(r,s,z)tZr,  (6) 

where  integration  is  over  the  whole  space.  Evaluation  of  the 
integral  in  Eq.  (6)  over  Eq.  (1)  leads  to 


where  UImn=p,s[Smn-p'mn/(2l+l)l  The  solution  Flmn(t)  has 
the  following  form: 

F'm [Blmno\i  exp(-  X‘t) ,  (11) 


/  =  1 ,2,3,4  for  Z  3=  2,  and  Z  =  1,2  for  /  =  1 , 2.  The  eigenvalues 
kj  for  l=?2  is  given  by 


-  (1/2)  j  (n00 + n22  ±  n0_0  ±  n2_2) 


(n'-n'2±n(; 


10-0 


n2-2) 


/  Re[n20] 
\Im[II20] 


(12) 


for  Z=l,2,  and  for  /  =  3 , 4,  the  sign  +  before  the  square 
brackets  in  Eq.  (12)  is  replaced  by  — .  For  /  =  0 , 1  two  eigen¬ 
values  are  X;=IIqo±IIq_0.  The  constant  coefficients  [R(„„  ],■ 
can  be  determined  using  linear  algebra  under  the  initial  con¬ 
dition  Eq.  (10)  [31]. 

In  SP  the  angular  distribution  component  of  Im 
=  [l .  QJJ .  V],  with  the  initial  polarized  state  Ispl0),  is  ob¬ 
tained  by 


PF(s,s0  ,t)/dt+  p,flF(s,s0,f) 


+  P-y 


F(s,s0,r)  -  J 


P(s,s')F(s',s0,t)t/s' 


I(0)^(s-So)<5(r-0). 


(7) 


The  gradient  term  disappears  in  Eq.  (7)  because  of  the 
Gauss-Stokes  law.  Hence,  F(s,s0,t)  can  be  exactly  solved  by 
expanding  in  GSF’s.  When  So  is  set  along  the  z  direction,  the 
initial  polarization  state  is  set  as  n0,  and  the  initial  reference 
plane  is  set  as  the  x-z  plane,  we  have 


Fmnn  (s,s0,f)  =  E 


FlmnSt)P‘  (COS  0) 


21+1  _  , 

-  Sm„—. - e  ft'Po0(cos  0) 


0  4  77 


Xexp(-  Zn0</))exp(-  p.at) , 


(8) 


Ff  (s,s0,f)  =  [T-‘F(s,  s0,  Z)TI^°)]m.  (13) 

The  angular  distribution  for  the  parallel  polarized  compo¬ 
nent  is  given  by  F||(s,s0,r)  =  [F2/>(s,s0,t)  +  Fg/)(s,s0,r)]/2. 
Similarly,  we  have  F±=(FSjP- -FsQp)/2,  FR={FSIP +FSVP)I2, 
and  Fl=(Fsip-Fsvp)/2. 


D.  The  cumulant  expansion 

The  polarized  photon  distribution  7m(r,s,f)  in  an  infinite 
uniform  medium  can  be  written  as  [31] 


4(r,s  ,t) 


where  (•  ■  ■)  means  the  ensemble  average  in  the  angular  space. 
We  perform  a  Fourier  transform  of  the  first  S  function,  then 
make  a  cumulant  expansion,  and  obtain 
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Im(r,s,t)=Fm(s,s0,t) 


(Sj5  /rfqexpjiq  .r 


(—  ic)k  ^r\ 

k=  1  k'  jk  h 


X 


I  dtk---  I 
Jn  J  o 


71^(4) 


./r 


(15) 


<R j)m„n  =  c[  dt'  fds'2 

Jo  J  n 

=  c2  |  ^L,„0-/cos  0)e-,<no-j)lf‘yj^  2 


Fm„(s,s',f-f')ijFrao(s',s0,?') 

2(1  -/?)+! 
47T 


^Dmn  nf(j)(l  -  h,\,n,0\l,n)(l  -  h.lpiQ,- j\l,nQ- j) 

21  +  1  0(cos  <9)<r'”<^  l  (19) 


- 


Y°  m’"0  4, 


where  T  denotes  time-ordered  multiplication  [38],  and  the 
subscript-index  c  denotes  cumulant.  The  cumulants  are  re¬ 
lated  to  the  moments.  For  an  arbitrary  random  variable  A,  we 
have  (A)C=(A),  (. A2)C=(A2)-(A)(A ),  and  so  on.  The  corre¬ 
sponding  moment  can  be  estimated  using  a  standard  time- 
dependent  Green’s-function  approach,  which  is  given  by 


- 1 

[  dtk~ 

j 

Jo 

Jo 

dt{Asjk(t^---su(tl)'] 


1 


Fm(s,s0,t) 


f  dtk  f  dtk_ ,  ■■  ■  f  dt  |  f  dslk) 

Jo  Jo  Jo  J 


X 


t-tyk) 


ik 


XF(SW,s^-1),4-  4-t)^.^  •  •  •  F(s(2),s(1),t2  -  h)sV 


with  y±1=  +i,  7o=l,  is  Clebsch-Gordan 

coefficients  in  angular  momentum  theory  [36],  presented  in 
the  Appendix,  and 


l’J 


exp(-  ht)  -  exp(-  \\t) 


K  -  \ 


l-h 


Xexp (~  flat). 


(20) 


where  i,j=  1 ,2,3,4  for  Z^2,  and  ij  =1,2  for  /=  1 ,2. 

In  SP  the  component  of  the  unnormalized  first  moment  of 
/„,  =  [/,  Q>  U,  V],  with  the  initial  polarized  state  Ixp(°l,  and  a 
=x,y,z,  is  obtained  by 


XF(s(1),s0,ti-0)I(0) 


+  (perm)  f , 


(16) 


<R XP  =  2  Uaj[T-l(Rj)TISpW]m.  (21) 

j 


where  the  abbreviation  “perm”  means  all  k\ - 1  terms  ob¬ 
tained  by  permutation  of  {/,},  i=l, ...  ,k  from  the  first  term, 
and  F(s6\s6-0;  ,)  is  the  exact  solution  given  by  Eq.  (8) 

in  Sec.  II  C.  If  we  cut  off  the  summation  in  Eq.  (15)  at  k 
=  2,  the  integration  leads  to  a  Gaussian  spatial  distribution 
characterized  by  the  first  and  second  cumulants. 


The  center  of  photon  distribution  for  the  linearly  polarized 
components  is  obtained  by 


KJ  M)  = 


[(Ra)f±(Ra)Sf] 

2Fp(s,s0,t) 


(22) 


1.  First  cumulant  (center  of  the  spatial  distribution) 

In  CP  it  is  convenient  to  set  s  in  a  spherical  harmonic 
basis. 


s  =  [hAoA-i]  =  [-  2  1/2  sin  0e+"f>, cos  0,  +2  1/2  sin  0e  '^], 

(17) 

which  is  related  to  Cartesian  component  basis  a=x,y,z  by 
sa=1,UajSj,  7=1, 0,-1  with 


where  the  +  sign  holds  for  the  p  =  \\  component,  and  the  — 
sign  holds  for  the  /?=_L  component.  The  center  of  photon 
distribution  for  the  circularly  polarized  components  is  ob¬ 
tained  by 


KJ  s’f): 


[(Ra)f±(Rg)f] 

2Fc(s,s0,t) 


(23) 


-  1 


U=  2~m  i 


0  1 
0  i 


0  2 1/2  0 


(18) 


Using  the  recurrence  relation  and  the  orthogonality  relation 
of  GSF,  the  unnormalized  first  moment  in  CP,  when  s0  is 
along  the  z  direction,  is  obtained; 


where  the  +  sign  holds  for  the  c=R  component,  and  the  — 
sign  holds  for  the  c=L  component. 

2.  Second  cumulant  (width  of  the  spatial  distribution) 

Similarly,  we  can  obtain  the  expression  for  the  unnormal¬ 
ized  second  moment  in  CP, 
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<r/2r/i>„M() = c22  \  p'm,n0-jrj2(c os  e)e-i(n°-ji  2  2  2 


2{l-h2-hl)+  1 
4t t 


n2  n  j  h2  hi 

X  -  h2,  l,n2,0|/,n2)<Z  -  h2,  l,n0  — 7i,— 72l^«o  “7i  -72><*  -  /;2  -  ^i>  l,«i.0|/  -  A2,"i> 


2/  +  1  r 


X(/-/i2-/71,l,H0,-;1|/-/!2,n0-7i)- -e  ^'p'o 0(cos  0)e 


(24) 


with 


i,j,f 


(  exp(-  x). ,i2  Alr)  -  exp(-  X-f) 

\  (X^2  -  a^-VAi)(x'  _  x}-*2-*1) 

exp(-  x(~*2f)  -  exp(-  X|f)  \ 

(. \‘rh 2  -  x}_,,2_/,1)(x'  -  x'-'12)  /  ’ 


where  i ,j ,/=  1 , 2 , 3 , 4  for  1 3*  2,  and  i  ,j  ,/=  1 , 2  for  /=  1 , 2. 

In  SP  the  component  of  the  unnormalized  second  moment 
of  Im  =  [I ,Q,U ,V\,  with  the  initial  polarized  state  is 

obtained  by 


<  =^2  \(uajuPh + t/„,2t/fti) 

A  22 

X[T_l(Rj2RJi)TI,s/,(0)]m.  (26) 

The  square  of  the  width  of  the  photon  distribution  for  the 
linearly  polarized  component  is  determined  by 


The  ballistic  component  has  been  subtracted  in  the  ex¬ 
pressions  for  the  angular  distribution  Eq.  (8),  for  the  first 
cumulant  Eq.  (19),  and  for  the  second  cumulant,  Eq.  (24). 
This  subtraction  ensures  that  the  summation  over  /  converges 
even  at  very  early  times,  as  discussed  in  Sec.  Ill  of  Ref.  [30]. 

In  the  case  of  the  transmission  of  light,  the  formalism 
produces  the  correct  first  and  second  cumulants,  but  the 
Gaussian-shaped  distribution  obtained  for  early  times  vio¬ 
lates  causality,  as  manifested  by  the  photons  at  the  front  edge 
of  the  distribution  traveling  faster  than  the  speed  of  light  in 
free  space.  We  have  introduced  an  approach  to  reshape  the 
distribution  in  Ref.  [30],  which  maintains  the  correct  center 
position  and  width  of  spread  of  the  distribution,  and  satisfies 
the  causality  condition.  This  technique  can  be  extended  to 
include  polarization,  combining  the  formula  for  the  calcula¬ 
tion  of  cumulants  of  polarized  components,  in  the  study  of 
transmission  of  linear  and  circular  polarized  light. 


III.  RESULT 
A.  Theoretical  result 


[(RaRg)r±(RjySf] 

2F,|(s,s0,f) 


-XA’  W 


where  the  +  sign  holds  for  the  p  =  \\  component,  and  the  — 
sign  holds  for  the  p  =  T  component.  The  square  of  the  width 
of  photon  distribution  for  the  circularly  polarized  compo¬ 
nents  is  obtained  by 


Dcafs(s,t) 


[(R„R;3)f±(RaRff)f] 

2Fc(s,s0,f) 


~K,X,b>  (28) 


where  the  +  sign  holds  for  the  c=R  component,  and  the  — 
sign  holds  for  the  c=L  component. 


F.  Gaussian  distribution  of  the  polarized  light 

Polarized  photon  intensity  7„(r,s,f),  v=  II ,  T  ,R,L,  is  ex¬ 
pressed  by  an  approximate  Gaussian  distribution. 


/„(r,s,f): 


Fv(s,s0,t) 


(4tt)3/2  [det  79  J 1/2 

X(ra-KJ(rf}-KJ  f  ■ 


exp]  2  -  7 [(°v)  '] 
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(29) 


We  now  use  the  formalism  developed  in  Sec.  II  to  study 
backscattering  of  circularly  polarized  light  in  a  scattering 
medium.  A  right  circularly  polarized  beam  of  light  with  ini¬ 
tial  Stokes  parameters  1=  [1,0, 0,1]  and  wavelength  X 
=  610  nm  is  assumed  to  be  injected  into  an  infinite  medium 
comprised  of  polystyrene  spheres  (with  refractive  index  m 
=  1.59)  suspended  in  water  (with  refractive  index  m0=1.33). 
The  absorption  coefficient  of  light  in  the  medium  is  assumed 
to  be  zero.  The  light  is  incident  along  the  z  axis  from  the 
origin  (x,y,z)  =  (0,0,0).  The  backscattered  light  is  received 
by  a  detector  located  at  (xrf,0,0).  Figure  1  shows  the  com¬ 
puted  time-resolved  profiles  of  circularly  polarized  backscat¬ 
tered  light  as  a  function  of  time  t  (in  a  unit  of  Ilrl  c),  collected 
at  direction  0=180°  and  position  xd=0.5lIn  where  l,r  is  the 
transport  mean  free  path,  for  polystyrene  spheres  of  diameter 
(a)  <7=0.10  /r m,  (b)  <7=0.213  /zm,  (c)  <7=0.855  /zm,  and  (d) 
<7=8.0  /zm.  A  remarkable  feature  is  that  the  backscattered 
signal  is  dominated  by  the  left  circular  polarized  (LCP)  light 
if  the  size  of  scatterers  is  small  compared  to  the  wavelength 
of  light,  as  shown  in  Figs.  1(a)  and  1(b),  but  it  is  dominated 
by  the  right  circular  polarized  (RCP)  light  if  the  size  of  scat¬ 
terers  is  large,  as  displayed  in  Figs.  1(c)  and  1(d).  Figure  2 
shows  how  the  peak  intensity  of  RCP  light  and  LCP  light 
vary  with  the  size  of  scatterers. 
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FIG.  1.  (Color  online)  Time-resolved  profiles  of  RCP  (solid 
curve)  and  LCP  (dotted  curve)  backscattered  light  intensity  as  func¬ 
tions  of  time  t  (in  the  unit  of  ltr/c),  received  at  angle  0=180°  and 
detector  position  xd=0.5Ilr,  where  l,r  is  the  transport  mean  free  path, 
for  different  diameters  of  scatterers:  (a)  0=0. 10  yum,  (b)  d 
=  0.213  yum,  (c)  0=0.855  yum,  and  (d)  0=8.0  /tm.  The  initial 
Stokes  parameters  I  =  [  1 , 0,0,1]  represent  a  RCP  light  beam  inci¬ 
dent  along  the  z  direction,  and  the  wavelength  \  =  610  nm. 


We  explain  this  phenomenon  by  considering  the  effect  of 
Mie  single  scattering  together  with  the  effect  of  multiple 
scattering.  The  Mie  formula  indicates  that  when  m/m0 
=  1.19,  for  small  particles  the  Stokes  V  component  is  nega¬ 
tive  for  0  larger  than  90°,  which  means  that  the  helicity  of 


FIG.  2.  (Color  online)  The  peak  intensity  of  the  time-resolved 
profile  of  RCP  and  LCP  backscattered  light  as  functions  of  0/X,  at 
angle  6=  180°  and  detector  position  xd=0.51tr.  The  incident  beam  is 
RCP. 


light  singly  backscattered  from  a  small  particle  will  flip.  As 
the  particle  size  increases,  the  helicity  of  the  near-backward 
scattered  light  computed  using  Mie  formula  becomes  oscil¬ 
lating,  and  the  region  of  negative  V  squeezes  to  a  range  near 
0=180°.  For  near-forward  scattering,  the  value  of  V  is  posi¬ 
tive,  hence,  the  helicity  does  not  flip.  For  large  particles  the 
anisotropy  factor  g  is  large,  and  the  probability  of  near¬ 
forward  scattering  is  high.  For  large  particles  compared  to  X, 
backscattering  results  from  the  accumulated  effect  of  many 
small-angle  (near-forward)  scattering  events,  each  of  which 
changes  the  direction  only  slightly.  These  small-angle  scat¬ 
tering  events  do  not  change  the  helicity  of  circularly  polar¬ 
ized  light.  Hence,  the  backscattered  light  is  dominated  by  the 
component  that  maintains  the  original  helicity.  When  par¬ 
ticles  are  small,  however,  the  events  with  large-angle  scatter¬ 
ing  play  a  more  important  role.  The  backscattered  light  is 
dominated  by  the  component  with  the  flipped  helicity. 

In  view  of  the  above  explanation,  it  is  expected  that  cir¬ 
cularly  polarized  light  with  reversed  helicity  dominates  only 
when  the  detector  is  located  close  to  the  source.  When  the 
distance  between  the  source  and  the  detector  becomes  large, 
more  backscattered  photons  come  from  the  accumulation  of 
many  small-angle  scattering  events,  and  those  retain  the 
original  helicity.  Figure  3  shows  the  computed  time-resolved 
profiles  at  different  source-detector  distances:  (a)  xd=0.1lln 
(b)  xd=0.2ltn  (c)  xd=l.0llr  for  a  fixed-particle  diameter,  d 
=  0.7  fx m,  and  detection  angle  0=180°.  We  see  that  the 
dominant  component  changes  from  LCP  to  RCP  as  the 
source-detector  distance  increases.  Figure  4  shows  the  back- 
scattered  time  profiles  for  different  detection  angles:  (a) 
cos  0=— 1,  (b)  cos  0=-O.9,  (c)  cos  0=-O.8,  for  particle  diam¬ 
eter  0=0.855  ytim  and  the  source-detector  distance  xd=  1.0 ltr 
While  total  intensity  increases  with  a  decrease  of  angle  0,  the 
RCP  component  remains  dominant,  which  indicates  that  the 
helicity  of  backscattered  light  is  not  sensitive  to  the  angle. 
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t  (units  of  xd  / c ) 

FIG.  3.  (Color  online)  Time-resolved  profiles  of  RCP  (solid 
curve)  and  LCP  (dotted  curve)  backscattered  light  as  functions  of 
time  t,  at  different  source-detector  distances:  (a)  xd=0.lllr,  (b)  xd 
=  0.2 ltr,  (c)  xd=l.0lln  for  particle  diameter  d= 0.7  fim,  the  detection 
angle  0=180°,  and  the  wavelength  \  =  610  nm.  The  incident  beam 
is  RCP. 

B.  Comparison  with  experiments 

We  compare  the  above  theoretical  result  with  that  of  our 
experiments  [23].  The  setup  for  the  time-resolved  back- 
scattering  experiment  is  schematically  shown  in  Fig.  5.  The 
sample  comprises  polystyrene  spheres  suspended  in  deion¬ 
ized  water  in  a  6  X  6  X  10  cm3  glass  cell.  Ultrashort  light 
pulses  of  100  fs  duration  are  generated  at  a  repetition  rate  of 
82  MHZ  with  wavelength  X  =  6 1 0  mn  by  a  colliding  pulse 
mode-locked  dye  laser.  A  quarter-wave  plate  is  used  to  obtain 
RCP  light.  Time-resolved  circularly  polarized  backscattered 
light  is  monitored  by  a  2  ps  resolution  streak  camera. 

Figure  6  shows  the  time-resolved  profiles  of  backscattered 
light  measured  separately,  with  (a)  small  particles 
d=  0.213  /rm  (g  =  0.389),  scattering  coefficients 
fislc= 0.61  cm-1  (/„.=  2.68  cm),  and  (b)  large  particles 
d=  8.0  ii m  (y  =  0.9 1 1 ),  scattering  coefficients 
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FIG.  4.  (Color  online)  Time-resolved  profiles  of  RCP  (solid 
curve)  and  LCP  (dotted  curve)  backscattered  light  as  functions  of 
time  t,  at  different  angles:  (a)  cos  B=  —  1,  (b)  cos  8=- 0.9,  (c)  cos  8 
=  -0.8,  for  the  particle’s  diameter  d= 0.855  /jl m,  the  detector  posi¬ 
tion  xd=l.0ltn  and  the  wavelength  X  =  610  nm.  The  incident  beam  is 
RCP. 

pis/c  =  0.61  cm”1  (/„.=  18.42  cm).  The  detected  backscattered 
light  is  dominated  by  the  LCP  (dotted  curve)  for  small  par¬ 
ticles  in  case  (a),  but  it  is  dominated  by  the  RCP  (solid  curve) 
for  large  particles  in  case  (b),  even  though  the  average  xd/llr 
for  case  (b)  is  much  smaller  than  that  for  case  (a).  Figure  7 
shows  a  comparison  of  normalized  circularly  polarized  inten¬ 
sity  of  theoretical  and  experimental  results  for  the  case  of 
<r/=  8.0  yum  (c//\=  13. 1  I).  Considering  the  differences  be¬ 
tween  the  theoretical  parameters  and  the  experimental  setup 
(not  an  infinite  medium,  uncertainty  in  source-detector  dis¬ 
tance)  theoretical  predictions  are  in  reasonable  agreement 
with  the  experimental  result. 
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FIG.  5.  Schematic  diagram  of  the  experimental  setup  for  the 
time-resolved  polarized  backscattering  measurement  (BS=beam 
splitter,  L  =  lens,  M  =  mirror). 

IV.  DISCUSSION 

A  detailed  theoretical  formalism  for  the  calculation  of  the 
distribution  of  polarized  light  intensity  scattered  from  an  in¬ 
finite  uniform  turbid  medium  has  been  presented.  The  for¬ 
malism  is  based  on  an  analytical  cumulant  solution  of  the 
vector  radiative  transfer  equation  for  an  arbitrary  phase  ma¬ 
trix.  It  is  then  specialized  for  a  phase  matrix  obtained  from 
the  Mie  theory  since  many  practical  applications  deal  with 
Mie- scattering  cases.  The  scattered  polarized  photon  distri- 


FIG.  6.  The  experimental  time-resolved  profiles  of  RCP  (solid 
curve)  and  LCP  (dotted  curve)  backscattered  light  as  functions  of 
time  t  for  (a)  small  particles  d=0.213  /rm  (g  =  0.389),  the  scattering 
coefficients  fJ,s/c= 0.61  cm-1  (llr= 2.68  cm),  (b)  large  particles  d 
=  8.0  fim  (g  =  0.911),  /jlsIc= 0.61  cm-1  (llr=  18.42  cm).  The  incident 
beam  is  RCP. 


FIG.  7.  (Color  online)  Comparison  of  theoretical  and  experi¬ 
mental  results  of  normalized  circular  polarized  light  intensity  for 
the  medium  with  a  particle  diameter  d=  8.0  fim  and  wavelength  A. 
=  610  nm.  The  incident  beam  is  RCP. 

bution  depends  on  the  size  of  the  scatterers,  the  distance 
between  the  light  source  and  the  detector,  as  well  as,  on  the 
detection  angle. 

The  advantage  of  the  analytical  cumulant  solution  is  that 
it  enables  the  fast  and  accurate  calculation  of  the  temporal 
profile  of  scattered  polarized  photon  distribution.  The  use  of 
Stoke’s  representation  allows  the  investigation  of  light  of  any 
polarization  state,  including  linear  polarization  and  circular 
polarization,  which  are  commonly  used  in  practical  applica¬ 
tions.  We  have  used  a  collimated  ultrashort  input  pulse  in  our 
formalism,  which  represents  the  experimental  conditions  us¬ 
ing  picosecond  and  femtosecond  lasers  more  closely  than  the 
plane  wave  cases  used  by  earlier  works  [1-6,8, 9].  In  our 
formalism,  the  linear  polarization  case  is  handled  by  setting 
the  initial  Stokes  vector  as  [1,  1,  0,  0].  Our  initial  results  for 
linear  polarization  are  in  good  agreement  with  earlier  results 
[4,5]. 

In  this  paper,  the  focus  has  been  on  the  calculation  of 
circularly  polarized  light  intensity  in  the  backscattering  ge¬ 
ometry.  The  formalism  enables  the  calculation  of  the  inten¬ 
sity  distribution  of  multiple  backscattered  light  components 
that  are  copolarized  and  cross  polarized  with  the  incident 
beam  for  different  angles  of  detection  and  for  different 
source-detector  positions.  It  is  applicable  for  the  analysis  of 
three-dimensional  spatially  resolved  measurements  of  char¬ 
acteristics  of  scattering  media,  such  as  satellite-based  lidar 
measurements  of  clouds  and  aerosols.  The  approaches  using 
a  plane-wave  incident  beam  calculate  the  backscattered  pho¬ 
ton  flux  at  the  plane  of  incidence  (,:  =  ()  plane),  and  simplify 
the  three-dimensional  probing  problem  as  a  one-dimensional 
problem  [2-6, 8, 9].  The  approaches  using  Monte  Carlo  simu¬ 
lation  can  also  handle  three-dimensional  spatially  resolved 
measurements,  but  require  a  long  computation  time 
[7,39,40], 

Once  the  copolarized  and  cross-polarized  scattered  inten¬ 
sity  distributions  are  calculated,  parameters,  such  as  the  de¬ 
gree  of  polarization  and  image  contrast  can  be  readily  calcu¬ 
lated  to  take  advantage  of  the  polarization  memory  effect 
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[25].  One  of  the  results  of  this  work  elucidates  the  depen¬ 
dence  of  the  relative  intensities  of  the  backscattered  copolar¬ 
ized  and  cross-polarized  light  on  the  size  of  the  scatterers  in 
the  turbid  medium.  We  find  that  when  the  incident  light  is 
right  circularly  polarized,  the  backscattered  signal  is  domi¬ 
nated  by  LCP  light  if  the  size  of  scatterers  is  smaller  than  the 
wavelength  of  light,  but  it  is  dominated  by  RCP  light  if  the 
scatterer  size  is  larger  than  the  wavelength  of  light.  The  con¬ 
clusion  is  in  good  agreement  with  our  experimental  results 
[23].  The  dependence  of  the  helicity  of  circularly  polarized 
backscattered  light  on  the  scatterer  size,  and  the  polarization 
memory  effect  have  been  investigated  by  other  authors 
[6,9,25-27,41],  Mackintosh  et  al.  [25]  experimentally  dem¬ 
onstrated  the  polarization  memory  effect  for  spherical  par¬ 
ticles  suspended  and  undergoing  Brownian  motion  in  a  liquid 
by  fitting  a  simple  diffusion-based  model  to  account  for  tem¬ 
poral  correlations  of  the  intensity  fluctuations  for  different 
polarization  channels.  The  theoretical  estimation  of  back- 
scattering  of  the  circular  polarized  light  using  the  numerical 
solution  of  the  VRTE  has  been  performed  by  Kim  and 
Moscoso  [6],  as  well  as  Sakami  and  Dogariu  [9],  which  has 
been  mentioned  earlier  in  the  Introduction.  Phillips  et  al.  [41] 
used  an  electric-field  Monte  Carlo  method  to  study  the  back- 
scattering  of  circularly  polarized  beams  normally  incident  on 
a  half  space  of  scattering  particles  and  found  that  the  back- 
scattered  light  of  the  same  helicity  formed  a  ring  centered  on 
the  point  of  incidence.  A  Monte  Carlo  approach  has  been 
used  by  other  researchers  as  well  [7,39,40].  Xu  and  Alfano 
[27]  used  a  random-walk  approach  to  analytically  estimate 
the  characteristic  depolarization  length,  which  is  obtained  as 
an  average  over  the  entire  detection  solid  angle.  All  these 
approaches  converge  on  the  general  result  that  the  circularly 
polarized  light  of  the  same  helicity  dominates  backscattering 
when  the  scatterer  size  is  larger  than  the  light  wavelength,  as 
does  the  formalism  presented  in  this  paper  and  experimental 
results  [23].  The  salient  feature  of  this  work  is  that  it  pro¬ 
vides  an  analytical  approach  for  the  calculation  of  polarized 
scattered  light  intensity,  an  experimentally  measurable  quan¬ 
tity,  as  a  function  of  time,  detection  angle,  and  source- 
detector  separation. 

The  dependence  of  the  polarized  backscattered  light  on 
the  size  of  the  scatterers  in  the  medium  has  many  potential 
applications,  including  biomedical  imaging  [13-15],  flow  cy¬ 
tometry  [18,19],  investigation  of  biological  cell  differentia¬ 
tion  [17,18],  subsurface  imaging  of  cracks  and  corrosion  be¬ 
low  paint  layers  [42,43],  lidar-based  remote  sensing  of 
investigation  of  cloud  and  aerosol  distribution  in  atmosphere 
[20,21],  and  the  imaging  of  targets  in  turbid  water  [22]. 

Biomedical  and  biophysical  applications  represent  an  in¬ 
creasingly  important  area  of  application  of  polarized  light 
scattering.  As  early  as  1976,  polarized  light  scattering  was 
identified  as  a  “new  biophysical  tool”  [6],  The  scatterers  in 
biological  materials  include  cells,  cell  nuclei,  mitochondria 
that  are  in  the  Mie  scattering  domain  (typical  size  varies 
from  5-20  /im  as  compared  to  the  wavelength  of  light  in  the 
visible  to  NIR  of  0.5- 1.5  /im).  Polarization  effects  in  scat¬ 
tered  light  have  been  used  to  study  bacterial  suspensions  in 
water  [16,19],  distinguish  between  a  number  of  leukocyte 
types  in  flow  cytometry  [44,18],  cell  differentiation  [17],  and 
polynucleosome  superstructures  [45].  It  has  recently  been 
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shown  experimentally  that  backscattered  light  may  be  ana¬ 
lyzed  to  obtain  information  about  the  size  distribution  of  the 
cell  nuclei  [15,46],  which  in  turn  could  be  useful  for  cancer 
detection  as  the  cell  nucleus  size,  shape,  and  distribution 
changes  with  cancer  progression.  In  biomedical  imaging  ap¬ 
plications  polarization  effects  in  scattering  have  been  shown 
to  be  useful  in  examining  skin  [47],  and  subsurface  struc¬ 
tures  in  the  prostate  and  other  tissues  [48]. 

Another  important  area  of  application  of  polarized  light 
scattering  is  the  remote  sensing  of  the  atmosphere  and  earth. 
Radar  technology,  the  workhorse  for  remote  sensing  of  the 
earth  and  the  atmosphere,  makes  use  of  radio  waves  and  its 
depolarization  to  obtain  information  about  remote  targets  and 
ground  flora  [49].  While  single-scattering  approximation  of 
radio  waves  is  reasonable  for  radar  applications  [1],  for  prob¬ 
ing  of  dense  cloud,  fog,  and  aerosols  in  the  atmosphere,  mul¬ 
tiple  scattering  of  light  needs  to  be  considered,  which  the 
formalism  developed  in  this  article  does.  Different  scatterers 
in  the  atmosphere  have  varied  size  distributions  [1,50],  such 
as  cloud  and  fog  droplets  (radii  generally  <100  /mi,  mean 
radii  typically  2.5-5  /mi),  aerosols  (typically  <1  /mi),  and 
hydrometeors  (typically  1  /im  or  higher).  Lidar-based  inves¬ 
tigation  of  cloud,  aerosols,  and  other  atmospheric  scatterers 
employ  visible  and  near-infrared  light  (typically  532  nm  and 
1064  nm)  [20,51].  Scatterer  size  to  wavelength  ratio  the  dic¬ 
tates  the  use  of  Mie  theory  for  the  study  of  light  scattering  by 
clouds  and  aerosols.  In  the  CALIPSO  (Cloud-Aerosol  Lidar 
and  Infrared  Pathfinder  Satellite  Observations)  mission  of 
NASA  [51]  one  of  the  tasks  involves  probing  with  linearly 
polarized  532  nm  and  1064  nm  nanosecond  duration  pulses 
and  the  analysis  of  the  temporal  profiles  of  backscattered 
light  for  assessing  the  vertical  distribution  of  clouds  and 
aerosols  in  the  atmosphere.  Pal  and  Carswell  [21]  have  ob¬ 
served  spatial  variations  in  the  polarization  properties  of 
multiple-scattered  light  backscattered  from  clouds.  Clearly 
these  types  of  experimental  data  and  those  obtained  by  ex¬ 
tending  measurements  to  even  shorter  pulses  and  circularly 
polarized  beams  are  fertile  grounds  for  the  application  of  the 
formalism  presented  in  this  paper. 
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APPENDIX:  CLEBSCH-GORDAN  COEFFICIENTS 
AND  rf-FUNCTION 

For  the  convenience  of  the  reader,  we  list  the  formula  of 
angular  momentum  we  used  in  this  paper. 

The  following  formula  for  Clebsch-Goldan  coefficients 
(l-h,  1  ,m,-j\ l,m—j)  are  useful  for  our  calculation. 
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where  the  row  index  (from  above)  /  =  1 ,0,-1  and  the  column 
index  (from  left)  /?  =  1 ,0,-1. 

For  obtaining  dl 11mn(x),  the  following  recurrence  relation  is 
used: 


.max(|m|,|« 


\x) 


1 

2max(|m|,|n|) 


[2  max(|/M 
(|/«  -  n|)!(|m  +  n\) ! 


1/2 


X  ( 1  -  x)  lm~"l/2(  1  +  x)  lm+"l/2 ,  (32) 


4„W-(/_i)[(/2_m2)(/2_n2)]l/2 

x  {(2/  -  l)[/(/  -  l)x  -  mn\d'~l (x) 

-  /{[(/  -  l)2  -  m2][(l  -  l)2  -  n2]}md!-2n(x)}, 

(31) 


for  l>max(\m\  ,|n|),  with 


and  dlmn(x)=0  for  /<max(|/w|  ,|n|).  For  m= 0,  /?  =  0  we  have 
<^0(x)=l,  dlQQ(x)=x,  and 


ldlo0{x)  =  (21  -  1 )xd'00l(x )  -  (/  -  l)c?o02(x) .  (33) 


[1]  A.  Ishimaru,  Wave  Propagation  and  Scattering  in  Random  Me¬ 
dia  (Institute  of  Electrical  and  Electronic  Engineers,  New 
York,  1997). 

[2]  R.  L.  T.  Cheung  and  A.  Ishimaru,  Appl.  Opt.  21.  3792  (1982). 

[3]  Q.  Ma,  A.  Ishimaru,  and  Y.  Kuga,  Radio  Sci.  25,  419  (1990). 

[4]  A.  Ishimaru,  S.  Jaruwatanadilok,  and  Y.  Kuga,  Appl.  Opt.  40, 
5495  (2001). 

[5]  A.  D.  Kim  and  M.  Moscoso,  Phys.  Rev.  E  64,  026612  (2001); 
SIAM  J.  Comput.  23,  2075  (2002). 

[6]  A.  D.  Kim  and  M.  Moscoso,  Opt.  Lett.  27,  1589  (2002). 

[7]  R.  Vaillon,  B.  T.  Wong,  and  M.  P.  Mengiic,  J.  Quant.  Spec- 
trosc.  Radiat.  Transf.  84,  383  (2004). 

[8]  Y.  Jiang,  Y.  L.  Yung,  S.  P.  Sander,  and  L.  D.  Travis,  J.  Quant. 
Spectrosc.  Radiat.  Transf.  84.  169  (2004). 

[9]  M.  Sakami  and  A.  Dogariu,  J.  Opt.  Soc.  Am.  A  23,  664 
(2006). 

[10]  A.  A.  Kokhanovsky,  Light  Scattering  Media  Optics:  Problems 
and  Solutions,  3rd  ed.  (Springer  and  Praxis  Publishing,  Chich¬ 
ester,  2004). 

[11]  D.  Bicout,  C.  Brosseau,  A.  S.  Martinez,  and  J.  M.  Scdhmitt, 
Phys.  Rev.  E  49,  1767  (1994). 

[12]  R.  Carminati,  R.  Elaloufi,  and  J.  J.  Greffet,  Phys.  Rev.  Lett. 
92,  213903  (2004). 

[13]  A.  H.  Hielscher,  J.  R.  Mourant,  and  I.  Bigio,  Appl.  Opt.  36, 
125  (1997). 

[14]  S.  G.  Demos  and  R.  R.  Alfano,  Appl.  Opt.  36,  150  (1997). 

[15]  V.  Backman,  R.  Gurjar,  K.  Badizadegan,  I.  Itzkan,  R.  R. 
Dasari,  L.  T.  Perelman,  and  M.  S.  Feld,  IEEE  J.  Sel.  Top. 
Quantum  Electron.  5,  1019  (1999). 

[16]  W.  S.  Bickel,  J.  F.  Davidson,  D.  R.  Huffman,  and  R.  Kilkson, 


Proc.  Natl.  Acad.  Sci.  U.S.A.  73,  486  (1976). 

[17]  W.  S.  Bickel  and  M.  E.  Stafford,  J.  Biol.  Phys.  9,  53  (1981). 

[18]  B.  G.  deGrooth,  L.  W.  M.  M.  Terstappen,  G.  J.  Puppels,  and  J. 
Greve,  Cytometry  8,  539  (1987). 

[19]  B.  V.  Bronk,  W.  P.  van  de  Merwe,  and  M.  Stanley,  Cytometry 
13,  155  (1992). 

[20]  D.  M.  Winker,  W.  H.  Hunt,  and  C.  A.  Hostetler,  Proceedings  of 
Laser  Radar  Techniques  for  Atmospheric  Sensing,  edited  by  U. 
N.  Singh  [Proceedings  of  SPIE  5575,  8  (2004)]. 

[21]  S.  R.  Pal  and  A.  I.  Carswell,  Appl.  Opt.  24,  3464  (1985). 

[22]  G.  W.  Kattawar  and  M.  J.  Rakovic,  Appl.  Opt.  38,  6431 
(1999). 

[23]  X.  Ni  and  R.  R.  Alfano,  Opt.  Lett.  29,  2773  (2004). 

[24]  G.  D.  Lewis,  D.  L.  Jordan,  and  P.  J.  Roberts,  Appl.  Opt.  38, 
3937  (1999). 

[25]  F.  C.  MacKintosh,  J.  X.  Zhu,  D.  J.  Pine,  and  D.  A.  Weitz,  Phys. 
Rev.  B  40,  9342  (1989). 

[26]  F.  C.  MacKintosh  and  S.  John,  Phys.  Rev.  B  40,  2383  (1989). 

[27]  M.  Xu  and  R.  R.  Alfano,  Phys.  Rev.  Lett.  95,  213901  (2005); 
Phys.  Rev.  E  72,  065601  (2005). 

[28]  W.  Cai.  M.  Lax,  and  R.  R.  Alfano,  Phys.  Rev.  E  61,  3871 

(2000). 

[29]  W.  Cai,  M.  Lax,  and  R.  R.  Alfano,  J.  Phys.  Chem.  B  104,  3996 

(2000). 

[30]  W.  Cai,  M.  Xu,  and  R.  R.  Alfano,  Phys.  Rev.  E  71,  041202 
(2005). 

[31]  W.  Cai,  M.  Lax,  and  R.  R.  Alfano,  Phys.  Rev.  E  63,  016606 

(2000). 

[32]  M.  Born  and  E.  Wolf,  Principles  of  Optics,  7th  ed.  (Cambridge 
University  Press,  Cambridge,  England,  2002). 


056605-10 


ANALYTICAL  CUMULANT  SOLUTION  OF  THE  VECTOR... 

[33]  S.  Chandrasekhar,  Radiative  Transfer  (Clarendon,  Oxford, 
1950). 

[34]  I.  Kuscer  and  M.  Ribaric,  Opt.  Acta  6,  42  (1959). 

[35]  J.  W.  Hovenier  and  C.  V.  M.  van  der  Mee,  Astron.  Astrophys. 
128,  1  (1983). 

[36]  B.  M.  Brink  and  G.  M.  Satchler,  Angular  Momentum  (Claren¬ 
don  Press,  Oxford,  1962). 

[37]  M.  I.  Mishchenko,  L.  D.  Travis,  and  A.  A.  Lacis,  Scattering , 
Absorption,  and  Emission  of  Light  by  Small  Particles  (Cam¬ 
bridge  University  Press,  Cambridge,  2002). 

[38]  R.  P.  Feyman,  Phys.  Rev.  84,  108  (1951). 

[39]  M.  J.  Rakovic,  G.  W.  Kattawar,  M.  Mehrubeoglu,  B.  D.  Cam¬ 
eron,  L.  V.  Wang,  S.  Rastegar,  and  G.  L.  Cote,  Appl.  Opt.  38, 
3399  (1999). 

[40]  S.  Bartel  and  A.  H.  Hielscher,  Appl.  Opt.  39,  1580  (2000). 

[41]  K.  G.  Phillips,  M.  Xu,  S.  K.  Gayen,  and  R.  R.  Alfano,  Opt. 
Express  13,  7954  (2005). 

[42]  J.  H.  Ali,  W.  B.  Wang,  P.  P.  Ho,  and  R.  R.  Alfano,  Opt.  Lett. 
25,  1303  (2000). 

[43]  J.  H.  Ali,  W.  B.  Wang,  R.  R.  Alfano,  and  M.  K.  Kassir,  Theor. 
Appl.  Fract.  Mech.  41,  1  (2000). 

[44]  G.  C.  Salzman,  S.  B.  Singham,  R.  G.  Johnston,  and  C.  F. 


PHYSICAL  REVIEW  E  74,  056605  (2006) 

Bohren,  in  Flow  Cytometry  and  Sorting,  2nd  ed.  (Wiley-Liss, 
New  York,  1990). 

[45]  S.  Zeitz,  A.  Belmont,  and  Nicolini,  Cell  Biophys.  5,  163 
(1983). 

[46]  L.  T.  Perelman,  V.  Backman,  M.  Wallace,  G.  Zonois,  R.  Mano- 
haran,  A.  Nusrat,  S.  Shields,  M.  Seiler,  T.  Hamano,  I.  Itzkan,  J. 
van  Dam,  J.  M.  Crawford,  and  M.  S.  Feld,  Phys.  Rev.  Lett.  80, 
627  (1998). 

[47]  R.  R.  Anderson,  Arch.  Dermatol.  127,  1000  (1991). 

[48]  W.  B.  Wang,  J.  H.  Ali,  R.  R.  Alfano,  J.  H.  Vitenson,  and  J.  M. 
Lombardo,  IEEE  J.  Sel.  Top.  Quantum  Electron.  9,  228 
(2003). 

[49]  M.  I.  Skolnik,  Introduction  to  Radar  Systems  (McGraw-Hill, 
New  York,  1962). 

[50]  V.  E.  Zuev,  A.  A.  Zemlyanov,  Y.  D.  Kopytin,  and  A.  V.  Kuz- 
ikovskii,  High-Power  Laser  Radiation  in  Atmospheric  Aero¬ 
sols  (D.  Reidel  Publishing,  Dordrecht,  1984). 

[51]  D.  M.  Winker,  J.  R.  Pelon,  and  M.  P.  McCormick,  Proceedings 
of  the  Conference  of  Lidar  Remote  Sensing  for  Industry  and 
Environment  Monitoring  III,  edited  by  U.  N.  Singh,  T.  Itabe, 
and  Z.  Liu  [Proceedings  of  SPIE  4893,  1  (2003)]. 


056605-11 


Determination  of  light  absorption,  scattering  and  anisotropy 
factor  of  a  highly  scattering  medium  using  backscattered 

circularly  polarized  light 

M.  Xu".  M.  Alrubaiee^,  S.  K.  Gayen*  and  R.  R.  Alfano* 

"Department  of  Physics,  Fairfield  University,  Connecticut,  CT  06824 
^Institute  for  Ultrafast  Spectroscopy  and  Lasers  and  Department  of  Physics, 

The  City  College  and  Graduate  Center  of  City  University  of  New  York,  New  York,  NY  10031 

Email:  mxu@mail.fairfield.edu 


ABSTRACT 

The  absorption  coefficient,  the  scattering  coefficient  and  the  anisotropy  factor  of  a  highly  scattering  medium  are  determined 
using  the  diffuse  reflectance  of  an  obliquely  incident  beam  of  circularly  polarized  light.  This  approach  determines  both  the 
anisotropy  factor  and  the  cutoff  size  parameter  for  the  fractal  continuous  scattering  medium  such  as  biological  tissue  and 
tissue  phantoms  from  depolarization  of  the  backscattered  light. 

Keywords:  absorption,  scattering,  anisotropy  factor,  backscattering,  depolarization,  circular  polarization,  tissue,  turbid 
medium 


1.  INTRODUCTION 

The  estimation  of  the  optical  properties  (absorption  coefficient  //  a  and  scattering  coefficient  fis)  and  the  anisotropy  factor 
( g )  of  a  turbid  medium  is  an  important  problem  in  atmospheric,  biomedical  and  hydrologic  optics.  Light  diffuses  in  a 
highly  scattering  medium,  whose  behavior  essentially  depends  on  the  reduced  scattering  coefficient  g,'a  =  / /,  s  ( 1  —  g).  The 
diffusion  of  light  remains  almost  the  same  when  varying  /is  and  g  yet  keeping  ///,  unchanged.  This  explains  why  g,a  and  /i's, 
instead  of  /ia,  fxs  and  g,  are  commonly  used  in  optical  imaging.  Light  absorption  and  scattering  coefficients  depend  on  the 
optical  properties  of  constituent  particles  and  their  concentration.  The  g-factor  measures  how  much  forward  peaked  is  the 
angular  distribution  of  light  when  it  is  scattered  once  by  a  particle  in  the  medium,  and  is  independent  of  the  concentration 
of  the  scatterers.  The  value  of  the  g-factor  correlates  strongly  with  the  size  of  the  scatterer.  Determination  of  the  g-factor 
then  provides  a  direct  probe  of  the  size  of  scatterer,  which  constitutes  one  useful  fingerprint  in  discriminating  normal  and 
diseased  tissues.1’2 

One  way  to  estimate  the  g-factor  is  to  use  an  iterative  method  that  consists  of  repeatedly  solving  the  corresponding 
direct  radiative  transfer  problem  by,  for  example,  Monte  Carlo  simulations,  with  updated  optical  properties  until  the  irra- 
diance  matches  the  measured  values.3  Explicit  analytical  equations  have  also  been  developed  (analytic  inverse  radiative 
transfer  method)  for  determining  the  optical  properties  in  the  simplest  problems. 4  In  both  approaches,  the  type  of  the  phase 
function  of  light  scattering  by  the  medium  is  usually  assumed  known  and  the  Henyey-Greenstein  phase  function  5  has  been 
commonly  used. 

In  this  paper,  we  report  on  a  scheme  of  detecting  light  absorption  coefficient,  scattering  coefficient  and  the  g-factor 
simultaneously  from  measurement  of  the  circularly  polarized  backscattered  light  at  different  positions  along  a  line.  Light 
is  obliquely  incident  on  the  surface  of  the  medium  with  an  incidence  angle  ~  70  degree.  Both  co-polarized  and  cross- 
polarized  backscattered  light  are  measured.  The  difference  of  the  backscattered  co-polarized  and  cross-polarized  light 
intensities  is  used  to  find  the  g-factor  of  the  medium.  Experimental  demonstrations  of  this  approach  to  obtain  the  three 
parameters  (absorption  coefficient,  scattering  coefficient  and  the  g-factor)  of  tissue  phantoms  are  presented. 
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Figure  1.  The  ratio  of  the  uncoiling  length  over  the  transport  mean  free  path  lx/lt  and  the  anisotropy  factor  g  of  soft  tissue  of  a  fractal 
dimension  4.5  and  of  varying  cutoff  size  parameter  X  =  2nnllnax/ X. 


2.  THEORETICAL  FORMALISM  AND  EXPERIMENTAL  METHOD 

Light  propagation  in  a  turbid  medium  can  be  described  by  radiative  transfer  theory.  Diffusion  approximation  to  radiative 
transfer  theory  is  a  valid  model  for  light  propagation  if  the  medium  is  highly  scattering  and  optically  thick.  For  a  polarized 
incident  beam,  light  propagation  in  the  medium  accompanies  its  depolarization.  The  characteristic  lengths  governing  light 
depolarization  in  turbid  media  have  been  obtained  earlier.  6- '  The  depolarization  of  circularly  polarized  light  is  found  to 
be  determined  by  the  uncoiling  length  l  x,  the  distance  over  which  the  degree  of  circular  polarization  of  the  beam  drops  by 
1/e.  The  uncoiling  length  for  a  turbid  medium  is  given  by  ' 

lx  =  l ujVln-i-  (1) 

where  Ax  =  max(A0,Ai),  A0  =  A.  ^  [S*S2]  dfi,  Ai  =  A.  sft  [S*S2]  nd/i,  7  =  f /^(ISil2  +  | S2\2)dn,  n  = 
cos  9,  and  Si^(9)  is  the  amplitude  scattering  function  with  9  the  scattering  angle. 8 

Light  scattering  by  biological  tissue  and  cells  can  be  described  by  the  fractal  continuous  random  medium  model.  1,2 
In  this  model,  tissue  light  scattering  is  determined  by  two  morphological  parameters:  the  fractal  dimension  D  f  and  the 
cutoff  correlation  length  lmax.  The  parameter  D  f  reveals  the  weight  of  small  vs  large  scattering  centers  in  the  sample  and 
Imax  measures  the  maximum  size  of  scattering  centers.  Light  depolarization  characteristic  lengths,  including  the  uncoiling 
length  lx,  in  biological  tissue  and  tissue  phantoms  have  been  investigated. 9  Fig.  (1)  displays  the  ratio  lx/lt.  and  the  g-factor 
versus  the  cutoff  size  parameter  X  =  2irnlniax/ X  for  a  fractal  continuous  random  medium  with  D  f  =  4.5  where  n  is 
the  refractive  index  of  the  medium  and  A  is  the  wavelength  of  the  incident  beam  in  vacuum.  This  figure  suggests  that  the 
anisotropy  factor  g  (and  /max)  can  be  determined  if  the  fractal  dimension  D  f  and  the  uncoiling  length  lx  can  be  found  in 
addition  to  the  transport  mean  free  path  lt. 


Proc.  of  SPIE  Vol.  6435  64350J-2 


Sample 


<  ► 

A./4  @800-nm 


Figure  2.  Schematic  diagram  of  the  experimental  setup.  M  stands  for  mirror,  L  lens,  and  /  is  the  focal  length  of  the  CCD  camera. 
The  A/4  wave  plates  and  polarizers  are  used  to  generate  a  circularly  polarized  beam  and  to  detect  co-  or  cross-polarized  light  at  the 
wavelength  of  800nm. 


Our  procedure  for  the  determination  of  /ia,  fis  and  g  is  as  follows.  First  //,,  and  ///  are  determined  from  the  diffuse  re¬ 
flectance  of  the  sample  for  an  obliquely  incident  beam  at  various  wavelengths  centering  at  A.  Second,  the  fractal  dimension 
Df  is  determined  from  /t' ( A)  oc  A:!  Dt,  Third,  the  uncoiling  length  lx  is  found  from  the  depolarized  diffuse  reflectance  of 
a  circularly  polarized  incident  beam.  And  finally,  the  anisotropy  factor  g  is  obtained  by  the  ratio  of  l  x/lt  (see  Fig.  1).  The 
first  two  steps  can  be  combined  to  retrieve  p's  and  Df  simultaneously  from  fitting  to  the  diffuse  reflectance  data  at  all 
wavelengths  together  assuming  the  power  law  for  //((A)  is  strictly  satisfied. 


For  an  obliquely  incident  pencil  beam  at  the  origin  on  the  surface  z  =  0  of  the  sample  and  propagating  in  the  direction 
(sin  a,  0,  —  cos  a)  where  a  is  the  angle  of  incidence,  the  initial  ballistic  motion  of  the  beam  within  the  medium  is  in  the 
direction  (sin  [3, 0,  —  cos  0)  where  n  sin  /3  =  sin  a.  The  diffuse  reflectance  at  the  position  (x,  y,  0)  on  the  surface  is  given 
by10 


R{x,y) 


lt  cos  f3 

47 xp\ 


(1  +  He&Pl)  exp(-/Lieff'Pl)  + 


lt  cos  /3  +  2 ze 
4t Tpl 


(1  +  PeSP'2)  exp (~pesp2), 


(2) 


where  peS  =  03pa/k,  Pi  =  \f(x-lt  sin /3)2  +  y2  +  0  cos2  /3,  and  p2  =  y/{x-  Tt  sin/3)2  +  y2  +  ( 2ze  +  Tt  cos 0)2. 

The  parameter  ze  is  the  extrapolation  length,  dependent  on  the  refractive  index  mismatch  at  the  interface. 11  Eq.  (2)  is  used 
to  fit  the  measured  diffuse  reflectance  along  a  line  in  the  x  direction  on  the  surface  to  obtain  p  a  and  lt. 

For  an  incident  circularly  polarized  beam,  light  depolarization  presents  one  extra  channel  that  light  leaves  the  co¬ 
polarized  state  to  the  cross-polarized  state.  Denote  the  co-polarized  and  cross-polarized  backscattered  light  intensity  for  an 
incident  circularly  polarized  light  as  I± ,  respectively.  The  depolarized  diffuse  reflectance  is  proportional  to  A I  =  /  +  L 
and  can  be  modelled  by  the  same  expression  (2)  with  the  only  modification  of  p,  a  — >  p!a  =  pa  +  l~l.  The  same  technique 
hence  can  be  used  to  fit  for  lx  in  the  region  far  away  from  the  incidence  point  where  the  diffusion  model  is  valid. 

The  experimental  setup  is  shown  in  Fig.  2.  The  light  source  is  a  Ti:sapphire  laser  (Spectra-Physics  Tsunami)  tunable 
over  the  750-860nm  spectral  range.  The  laser  was  operated  in  CW  mode  by  turning  off  the  acousto-optic  modulator  and 
had  a  spectral  FWHM  (full  width  at  half  maximum)  bandwidth  less  than  lnm.  The  beam  passes  through  a  converging  lens 
(focal  length  40cm)  and  focused  on  the  surface  of  the  sample  after  elevated  using  a  periscope  mirror  configuration.  The 
incidence  angle  is  set  to  be  a  =  72.4°.  The  refractive  index  of  the  diluted  Intralipid- 10%  suspension  is  n  =  1.334  and 
the  extrapolation  length  is  computed  to  be  ze  =  1 .69/(.  ' 1  Backscattered  light  in  the  direction  normal  to  the  surface  of  the 
sample  is  recorded  by  a  CCD  camera  for  the  wavelengths  from  760nm  to  860nm  with  a  step  of  lOnm. 
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Figure  3.  The  wavelength  dependence  of  the  absorption  and  the  reduced  scattering  coefficients  of  the  Intralipid- 10%  suspension. 


3.  RESULT 

Fig.  3  displays  the  wavelength  dependence  of  pa  and  //.(  from  fitting  the  total  diffuse  reflectance  where  light  is  collected 
regardless  of  its  polarization.  The  absorption  spectrum  follows  essentially  that  of  water.  The  wavelength  dependence  of  the 
reduced  scattering  coefficient  follows  p(,  oc  A-15.  The  fractal  dimension  of  the  Intralipid- 10%  suspension  is  then  4.5.  It 
should  be  noted  this  experimental  value  of  D  f  is  a  bit  larger  than  the  value  (Df=  4.2)  obtained2  from  the  Mie  scattering 
computation  based  on  the  particle  size  distribution  of  the  Intralipid- 10%  suspension  measured  by  Staveren  et.  al.  12 

The  depolarized  diffuse  reflectance  is  then  measured  at  the  wavelength  of  800nm  at  which  the  absorption  coefficient 
and  the  transport  mean  free  path  have  been  found  to  be  pQ  =  0.003mm_1and  lt  =  1.77mm,  respectively,  from  Fig.  (3). 
Fig.  4  displays  the  total  backscattered  light  /  =  /++/_  and  the  depolarized  backscattered  light  A I  =  I+  —  Nearby 
the  origin,  the  intensity  of  the  cross-polarized  light  is  stronger  than  that  of  the  co-polarized  light  and  A I  <  0.  In  the  far 
zone  (the  region  far  away  from  the  origin),  the  intensity  of  the  co-polarized  light  is  stronger  than  that  of  the  cross-polarized 
light  and  A I  >  0.  Fitting  of  the  depolarized  diffuse  reflectance  ex  A  /  along  a  horizontal  line  3 1 t  away  from  the  incidence 
point  inside  the  far  zone  yields  the  effective  absorption  fi'a  =  0.47mm  .  The  uncoiling  length  is  then  2.13mm  and  the 

ratio  lx/lt  =  1.20. 

From  the  analysis  of  the  depolarization  property  of  the  fractal  continuous  medium  of  D  f  =  4.5  (see  Fig.  1),  both 
the  g-factor  and  the  cut-off  size  parameter  are  now  obtained.  We  find  g  =  0.68  at  the  cutoff  size  parameter  3.4  with 
lmax  =  0.32p.m.  The  values  of  the  g-factor  and  the  cutoff  size  of  scatterers  were  reported  to  be  0.64  at  800nm  and  0.34pm 
in  Staveren  et.  al. 12  Their  agreement  demonstrates  the  effectiveness  of  our  proposed  method. 

4.  DISCUSSION  AND  CONCLUSION 

For  a  range  of  values  of  lx/lt  in  Fig.  1,  there  may  be  two  sets  of  the  g-factor  and  the  cutoff  size  parameter  giving  the  same 
ratio  lx/lt-  It  does  not  pose  any  difficulty.  The  left  branch  of  the  l  X/U  curve  when  the  cutoff  size  parameter  curve  is  less 
than  2.6  corresponds  to  I+  <  /_  and  its  right  branch  corresponds  to  I+  >  /_  far  away  from  the  incidence  point. 7  The 
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Figure  4.  The  intensities  of  backscattered  light  (a)  the  total  intensity  J+  +  /_  and  (b)  the  depolarized  intensity  /+  —  /_.  The  circularly 
polarized  light  is  incident  on  the  origin  and  with  an  incident  angle  of  72.4°.  The  CCD  camera  measures  I±  directly. 

experiment  shows  that  I+  >  /_  in  the  far  zone  (see  Fig.  4b)  and  the  right  branch  of  the  lx/lt  curve  should  be  used  in  our 
data  analysis. 

The  absorption  coefficient,  the  scattering  coefficient  and  the  anisotropy  factor  of  a  highly  scattering  medium  are  deter¬ 
mined  using  the  diffuse  reflectance  of  an  obliquely  incident  circularly  polarized  light.  This  approach  obtains  simply  both 
the  anisotropy  factor  and  the  cutoff  size  parameter  for  the  fractal  continuous  scattering  medium  from  circular  depolarization 
of  the  backscattered  light. 
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